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NOTE ON THE ALGEBRA OF ELLIPTIC FUNCTIONS 


By SIR CHARLES DARWIN (Newnham Grange, Cambridge) 


[Received 27 August 1959] 


SUMMARY 

Much time can be wasted over the algebra of elliptic integrals, and the present 
note gives a method which has proved very convenient in shortening this work. It 
has the advantage that every result can be derived by itself very nearly from first 
principles. 

It treats of the Jacobian functions, making use of their associated # functions, 
by transforming the products of pairs of such functions. 

The method can be widely extended to much more complicated cases, but here 
only familiar results are given, including the quadratic relations for the Jacobian 
functions, their derivatives, the addition theorems, and the relation #; = #, 3; A. 

The main formulae are given for elliptic integrals of the second kind. 

Examples of Landen transformations are exhibited. 


1. CoMPUTATIONS associated with elliptic integrals are apt to be intricate, 
and any simplified method is useful.t It is not claimed that the method 
here described is new, but it seems to be less known than it deserves. 
Here I shall only demonstrate some of the more familiar theorems, but 
experience shows that the method can be very widely used. It has one 
merit that is not found in some of the developments of the theory, which 
is that each proposition can be deduced by itself from first principles, 
instead of coming at the end of a long sequence. A good many propositions 
will be demonstrated, and to save space I shall sometimes refer back to 
results already proved, but it will be observed that each section of the 
paper could stand by itself, in that each of the relations it involves can 
be quickly derived just when it is needed and not in advance. 

There is a rivalry between the Weierstrassian and the Jacobian functions. 
Of these Weierstrass’s 9 has a certain formal elegance for complex func- 
tion theory, but there is no doubt that for practical purposes it is very 
inferior. Even if the first stage of its use may look simple, it will usually 
be found that this is later paid for by much trouble in evaluating the 
two complex periods. With the Jacobian functions this trouble is avoided 
at the start because one period is real and the other pure imaginary, and 


+ I assume some knowledge of the functions; for example the broad elements of what 
is given in Whittaker and Watson’s Modern Analysis, chs. xxi and xxii, but without 
attention to the details of proof. The present method excludes infinite products from con- 
sideration. 

(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 2, 1960] 
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for practical purposes in evaluating integrals, it is real values that are of 
main interest. 

Here then only the Jacobian functions are considered, and they are 
discussed with the use of the #-functions. In fact the processes to be 
used are equivalent to Landen transformations, but whereas the ordinary 
presentation of those transformations involves rather tiresome conversions 
of the variables of integration, here they come out as matters of elementary 
algebra, and indeed offer themselves as the natural steps to be taken in 
that algebra. 


2. The four 3-functions are 
B(x) = i > ( ugiu ie in(2u +1) 
Soar) = F qurPeineusn, 
ota dans ai (2.3) 
hy (x) > (—)uqieize™, (2.4) 


x 
In the present work > will always mean > , and, since q < 1, the con- 
x 


vergence of the q-series is always so strong that it needs no consideration 
at all. 
The relations of these functions to the elliptic functions are given by 


setting : ‘ 
é= rh 


and then 


3, 3, (2x) 
i, 94(a)’ 
3, 9,(2) 
Fy Pg(x)’ 
F, F(x 
dné = — al) 
Fs F4(x) 
Here the modulus k = 93/82 and hk’ = 34/83. 
There will also occur a similar set of 3-functions, but with argument 
q? instead of g. These will be denoted by using ¢ instead of 3. Thus, 
for example, 


ané = 


ené = 


d(x) = > (—)"q?"e2", ete. (2.9) 
In many of the cases given below the series vanish because their terms 
cancel in pairs. The reader will easily get practice in seeing these results, 
but for example the following all vanish: 
> (—)eqgu™,  S (2u4+- 1g", > (—)“ug’’, > ug’. 
To help in following the work I shall sometimes indicate these vanishing 
terms by inserting a 0 in the next expression deduced. 
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It is also very easy to check against blunders by expanding the q-series 
to the first two terms or so. 


3. The principle of the method is to take products of pairs of #-func- 
tions and to transform them into products of pairs of ¢-functions. 

The products must conform to one of two rules, which I shall call 
the [34] and the [12] rules. For the [34] rule the product must involve 
q’*", which implies that it is of one of the forms 3,( )3,( ), or A( )I4( ) 
or 3( )d,( ). This product is transformed by the substitution 


pS onde TS hes a 


which gives in the first term 
u?+v? = 2m?4+-2n? 
and in the second term 
u2+ v2 = 2(m+4)*?4+-2(n+ 4). 
In the [12] case only 3,( ) and #,( ) may occur, and a different substi- 
tution is more convenient. This is: 


=m-n u“ m+n 
Ze >. (" mea) ys es al (B) 


Then, in both terms here, 
(w+-4)?+(v+4)? = 2(m+4)?+2n?. 

As will appear, the results always give greene of series in which either 
¢, and ¢, go together, or else ¢, and ¢,. From this fact it is possible to 
work out results for the [23] and the [14] pairs by a ‘reverse’ transforma- 
tion, in which the work starts with the ¢’s instead of the #’s. The process 
is only a little more troublesome, and examples are given below. 

The third types of product [13] and [24] must not be used. 


yg 


4. A first ~— example illustrates the process. Using (A) 
# : > (- yuerg ut+e? _ Dd S [qrur* +n? — grin td? +0 +47) 
= 8-4 


Similarly, w= 62+ 43. 


Again, using (B), 
ie %y +4" +4 + 26 
# _ > ¢" yPre(rs4yP — +> qum+? 2n 2. 
since both terms give the same result, so 
2_« 
H = 224s. 
From these results it follows at once that 
4) 94 _ ga 
I+ df = I. 
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It will be noted that 
#3 9-8} _ 1 
a oo, ae 
which is the modulus corresponding to the Landen transformation. 
5. The important result 
H = 0,8,9, (5.1) 
is a deeper proposition than most others of elliptic function theory, but 
with the present method it is just as easily proved as the rest. Thus, 
applying (B), 
oe ty = , > 
¥. 


od 


( y“(2u j 1 gu tb? He +b? 
(- 


m+" (2m + 2n+ 1)gXm +d? +2n*2 
since both terms are identical, and this is equal to ' 
2(¢1 644-9). (5.2) 
\gain 
0,8, = > 2 (—)*g = BD [(— reg ts (— rte ign ihe) 
65 +0. (5.3) 
Then with the use of (4.3): 
oe _ 241 4 


— 1M _ f(g*) = fiq*) = fig? 
I, 39, , 2.3.63 f(q?) f(q') f(q ) 


f(q) 


Pt 
lim f(q) lim = 1 


q—0 q-0 =q*. bel 
which establishes the result. 
6. The results are just as easy to get when 2 does not vanish. Thus 
a, d,(a) > > ( yu ; tq ' v, ix2u 
a ig ¢ 2n*y ix(2m +4 2n)__g2im +4)242(n +h)*¢ ix(2m42n4 >| 
—_ 
3(x)—d2(x), 
using (A). Similarly, 
9 ' 9 
t, F5(x) b5(x) + 3(x). 


Again using (B) 
8, 8,(z) = > > qiutiProriletasu+y) 


DV OS g2im+4? 4 2n?pix(tm4+2n41)9 
ba han Tl 


== 26,(2)¢,(2). 
Thus at once 9% I2(x) +93 (x) = 92 93(zx). 
Dividing by 3j(x) and writing = x8, this gives 


, 2? 5 a9 
4 935 cn*é = 93 gine 
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or dn?é = k’?+ ken, (6.5) 
ay , r 
=z, and k’? — J 
3 v3 
The corresponding relation for sné is most easily derived by writing 
«+4n for x in (6.4), which yields 


92 92x) +93 92x) — 92 9%(z). (6.7) 


verifying that (6.6) 


To obtain the same result by the present method directly requires the 
evaluation of #, 3,(x), ete., which is also quite easy. 
An example ofa more complicated relation is 


9,(2)9g(y) = ESE (—)ngln edhe edtetadus sinaesy 
= bs > ( 7 ym+ngrm +4)? ron" gfatten +2n+1)+iy(2m—2n+1)_,_ 
| gix(dm+2n +1) +iy2n—2m Wy 
<5 TD (— ym engtion 8 204 piles yeam 41) +ir— em 4. 
{gil sian +ir—uXem +1] 
hy (e+ y)by(r—y) + by(2 +-y)b,(x—y). (6.8) 
7. The differentiations of the elliptic functions themselves are easily 
treated. From the [34] rule the best one to attack is dné. Thus 
F(x) g(x) — F(x) 9,(2) = YY (—)ugh’ etr2u+iz2e] (2y—i2u] 
~~ > > [( _)m +Mg2m* +2n*pirdm)y | 
{.(—)m+n +g im +4)? +2(n +4) ptridm *2(2n ' 1)] 
2i[ h4(2x).0— id, (22)d}] 
2; $1(22). (7.1) 
Also from (6.8) FB (x)I,(2) = b44,(22). 
So BH (x)F,(x)—P(x)94(2) = —2bob564¢,(22) = —BZH,(x)I,(x). (7.2) 
Then 


Ld dy F(x) | 


Is go I, (x)F,(x) , 
: -- — Da ar 3s — . 2 : B 
H dx I, H,(x) HR” F(x) ksn€eng. (7.3) 


né = 


The other differential relations are easily obtained by differentiating 
(6.4) or (6.7), but I will illustrate the ‘reverse’ process by finding the 
derivative of sné directly. This requires the evaluation of 

(1) hq(x)—h4(x)P (x) = J D (—)v regu td? serena sd +2 Iu + 1 — Qo). 
If (A) is applied to this the first exponent of q is 


4m? +-4n?-+-2m+-2n+-4 = (2m+4)?+(2n+ 4)? 
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and the second is (2m + $)?-+-(2n+ $)?. 
Thus the expression becomes 
_ / [qe +4 (2n +4)" piz(4m + V(4n a | ) qm +1? 12n+9"piridm +3)(4n P 3)! 

Now the expression 4{(—-1)"+-(—1)"] is +-1 if u and v are both even, and 

| if both are odd, and it vanishes if one is even the other odd. 

So 
f;(x)b,(x) $4(x) )y (x) >») 4 ( : 1)"+ "| +h 4404 einlusy +1)(2 v- 1) 

404d, (xh | bie I, 9,(x) = djdbo(x)b,(x) (7.4) 

by the use of (6.3). Since q is quite arbitrary, this equation in terms of q? 
will be equally true as an equation in q. So the ¢’s can be all changed into 
3's, and it easily vields the result 


wane enédné. (7.5) 


8. The addition formulae are a little more troublesome, but mainly 
because they can be cast in a great variety of forms only a few of which 
are familiar. 

The simplest procedure is as follows: 
from (6.8) F(x)F(y)+-F,(x)F,(y) = 26,(2+-y)dby(a—y). (8.1) 
Again, with a slight change of method so as to show the flexibility of the 
procedure, 

Ha(ar)F,(y) j I,(xr)Fg(y) > > qh +r? pirdus wan ( | yu | ( 1)"}. 
Applying (A) the last factor makes the sum of the second series vanish, 
and the first series is doubled. The result is therefore 

Y)by(a—y). (8.2) 
Then 


I (x)I(y) + F(x) Cy) b(e+y) _ 3 (r+ y)t, (8.3) 
s(x )4(y) + F4(x)Fa(y) d,(x+-y) “7, Pala } yd, + +, (2-4 YP 


Changing over into elliptic functions, and adopting an obvious notation, 
this gives k(8,C, +8 Cy) kS (8.4) 
d,+d, ~ ep’ ; 

which is one form of the addition theorem. 
To derive the standard form, observe that the right-hand side is 


(1—D)/(kS), and it is easy to change round the left side similarly. Then 


i—D (8; C2 +82 ¢4)(dy- -d;) 
ro @ (8; C5 -8_C det-d,)’ 
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So D—*! C2d,—82C, dy 
8, Cgd,—80¢,d, 
and this is easily reduced to the standard form. It is also easy to eliminate 


D and get the standard form for S. 
9. Elliptic integrals of the second kind can also be treated in the same 
manner. First there are curiously simple relations that are useful. 


aw 7” 


l i 2 2) 
5,8, ee te 
4 3 3”“4 
os -% > > [(—)m+mq2m*+2n*2m | 2n+ 
d 
+(— Jit +1 2m +4)? +2n +4 2am j 1)(2n +1)] 


04,2] = (2449)? = Of. (9.1) 


BR 


By differentiating (6.4) twice and putting 2 = 0, it is easy to see that 


consequently RNG 
—_ =H 9.2 
d, 4 ( ) 


and rr ~3 = Hf. (9.3) 
a Ys 

Alternatively the ‘reverse’ process easily proves (9.2) by working out 

$;¢.—¢5¢3. In this case the ‘reversing factor’ differs from that in section 

7, being here 3{1-+-(—1)"+?}. 


Again 9" 


98 2 geo] — ut 1) 
U's 2 


1 > > 2 2, an3 9 
EE et cere qm +4? +201 2m + 2n + 122 
2b2 3 


— | 14° 4.40442 4,] = 22499 — 2% 
= 5 4, eabs 0+¢3 $2] = at tg 


10. For the integral of the second kind evaluate 
F5(x)I (x) —F P(x) = ¥ DY (—)e req’ Petrrurizeey — 4y2+ suv] 


+44 


or better 
= > D(—)sttqy'trefenusey — 2)(u®+ v?—2ur) 
we EE an eran + i) 
= 265 4(2)—243 $422). pied 
Again F(x) = $3 $3(2x)—b, $4(22), 
3(2) = p35 ¢5(27) +4 ¢,(22). 
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which is left to the reader to verify. Then 

d F(x) - 1 (45 

dx d,(x) #(x)\d, 

= gp 282) (te #) 

B(x) © \dbs | be 


« 


[93() +.93(x)] — £3 9%(x)—93%(x)]| 
2 J 


, wv 
#4 dn2é +- —*. 
- ‘> 


Ve 


Integrate this with respect to 2 from 0 to 47, and it gives 


- o. 
0= 3 dn?é dx+ jn 


2 
(10.3) 


Integrating instead from 0 to x, 


r 
. 


» 
a | dnt dx — = 9B .x 
7 


F(x) ” 
By(x) q 
0 


» 
OE (¢) —= BE 


so that Z(é) = E(é)— BY oo 942) : 
BF O(a) 


Finally fi elliptic integrals of the third kind there is needed 
B b 

oO (x 
Zé) dé = | 28, 

J (x) 


x 


ix, where a = a2, B = bd} 


3 ,(b) 
#,(a) 


which is the well-known way of handling these integrals. 


11. When it is required to transform a cubic or quartic square-root into 
a form expressible by standard elliptic integrals, there is a prescribed 
routine for doing so. It consists in finding the two real factors, either one 
or both of them quadratic, and combining a sum of them so as to make 
a perfect square. This provides the transformation which makes the 
reduction always possible, but it will often be found not to be the best 
or simplest substitution. After carrying it out it is always worth testing 
whether a Landen transformation one way or the other may not yield 
simpler formulae. 
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Though all the work of this paper has been implicitly dealing with 
Landen transformations, explicit examples may be given in conclusion. 
é}_ 1k’ 
$2 14k” 


If ] is the modulus of the transformed functions, then / = 
. o, 1 6. a] 2B) ’ 

2 E ° — Pine es — m — grey 4 = i—f l), 

@) BO) = —deges = — raed git tl] = Ter OKO 

(3) dng = 26Fal*) _ Sz) — Pilz) _ 1—Ten'' 

"Pg Ay(x) F(x) +-GF(x) 1+ I sn’2E" 
where sn’ is the function to modulus /, and 
—— rd5 a g 
iil oa iif 84d 

From this the functions sné and en€ can also easily be expressed in 

rational forms in terms of sn’, en’, dn’. 
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BY PERTURBATION METHODS 
By P. LANCASTER 


(Department of Mathematics, University of Malaya, Singapore) 
| Received 24 September 1959] 


SUMMARY 

The problem of small vibrations of a conservative mechanical system is formu- 
lated. A method is developed by means of perturbation theory for the estimation 
of frequencies, rates of decay, and modes of vibration of the system when small 
damping forees are introduced. The analysis is based on a prior knowledge of the 
natural frequencies and principal modes of the undamped system. Special attention 
is given to the cases in which some natural frequencies of the undamped system are 
exactly equal, and those in which oscillations occur about a position of equilibrium 
which is critically stable. When several of the natural frequencies are nearly equal, 
a method of solution is given which can be applied when the variation in these 
frequencies is related to the magnitude of the damping forces in a particular way. 


1. Introduction 

THe problem posed by the small vibrations of a conservative holonomic 
system with a finite number of degrees of freedom is so well known that 
it requires very little introduction. The work of Lord Rayleigh (1) is still 
a standard reference work on this subject, and provides the basic theory 
for the developments of this paper. A recent paper by Bishop and Johnson 
(2) reviews this theory and derives some approximate results giving the 
effect of light viscous damping on the modes and frequencies of vibration 
of the undamped system. In this paper we considerably extend the results 
of Lord Rayleigh and of Bishop and Johnson. This extension is facilitated 
by employing as far as possible the powerful and suggestive methods of 
matrix algebra. 


In applying the perturbation method to any problem it must be as- 


sumed that all the properties of the unperturbed system are known. It is 
therefore assumed throughout this paper that all the natural frequencies 
of the undamped system are known, together with the corresponding 
principal coordinates. This is a well-known problem and has been 
thoroughly investigated by several authors. It is obvious that these 
properties of the undamped system must be known to considerably greater 
accuracy than the magnitude of the perturbations to be permitted. 

The perturbation method has been applied to matrix latent root 
problems by Jahn (3) in improving approximate latent roots and vectors, 


(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 2, 1960) 
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and by Van de Vooren (4) in investigating the effects of small changes 
in a parameter on the flutter characteristics of an aircraft. The prob- 
lem in which nearly equal latent roots (or natural frequencies) arise, 
is approached via that of exactly equal latent roots. The latter problem 
is solved by a method similar to that used by Courant and Hilbert (5) in 
applying perturbation theory to the analysis of self-adjoint linear differen- 
tial equations. 


2. The undamped system 

We consider a conservative holonomic system performing small oscilla- 
tions about a position of stable equilibrium with no external forces acting 
on it. In this case the kinetic and potential energies are positive definite 
homogeneous quadratic forms of the generalized velocities, ¢,, and coordi- 
nates, g,, respectively; where / = 1, 2,..., n and n is the number of degrees 
of freedom of the system. As in (1) and (6), the equations of motion can 
then be written in matrix forms, 


Aq+Cq = 0, (1) 
where A and (' are real symmetric matrices of constants, and q is the 
column vector with elements q,, qo,..., 7,- If we look for solutions of the 
form q = q,e”, we obtain 

(2A +C)q = 0, (2) 
(where we drop the suffix 0). The condition that these homogeneous 
equations in the q, should be consistent is |A24+-C) = 0, which shows 
that there are n possible values for A? such that 
(A7A+C)q, = 0, fori = 1, 2,..., n. 


This relation defines the latent roots A? exactly, but the latent vectors 
q; are not defined to within an arbitrary scalar multiplier. If we form the 
(n * n) matrix Y, whose ith column is the vector q;, and the (n * n) diagonal 
matrix A*, whose ith diagonal element is A?, then we have 


AQA?+CQ = 0. (3) 


It is shown in Appendix I that if the matrix Q is non-singular we may 
premultiply equation (3) by the transposed matrix, Q’, and obtain 


AAM+0E = 0, (4) 


where A = Q’AQ and C = Q’CQ are diagonal matrices. The above condi- 
tion on Q implies that the n latent vectors q; (which define the modes of 
vibration) should be linearly independent. It can be shown that this is 
always true in this case provided that either of the matrices A, C has 
rank n. Now this is in fact the case for matrix A (because it is derived 
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from the kinetic energy which cannot be negative). Thus we may permit 
C to have rank less than n and the diagonalization process may still be 
performed. If this is so, then some of the roots A; will be zero. This 
corresponds physically to oscillations about a position of equilibrium 
which is neutrally stable: i.e. in a mode of vibration corresponding to a 
latent root which is zero, the displacement may be proportional to the 
time. 

It can also be proved that all the non-zero latent roots of equation (2) 
are real and negative. Thus, if A? w?, the solutions of equation (1) 


will be of the form _ 
G; = do;(cos w, t+i sin w,t), 


where qo; is real. The real and imaginary parts of this expression also 
satisfy equation (1), and these give the solution physical significance. 


3. Natural frequencies distinct 

For the time being we will assume that the latent roots, A? of equation 
(2) are all distinct and non-zero. 

When small viscous damping forces are present in the system equation 
(1) must be replaced by 


Ag + 6q+Cq = 9, (5) 
where 8 is an (n <n) symmetric matrix whose elements are small com- 
pared with those of A and C. Thus, we view the entry of the term fq as ¢ 
small perturbation of the system discussed in section 2. For convenience 
we now make a transformation to principal coordinates of the undamped 
system, Py, Po... P,, given by q = Yp; and further, we suppose that the 
columns of matrix Q are normalized in such a way that A = J. Equation 
(4) then gives C \2. We use A?; for the diagonal elements of A? and 
reserve A; and A for the latent roots appropriate to the damped system. 
If we premultiply equation (5) by Q’ after making the above transforma- 
tion we obtain p LAp—Azp =k (6) 
where 8 = Q’8Q, and is clearly a symmetric matrix. We now write 
B = ¢B, say, where « is the perturbation parameter, and the elements of 
matrix B are of the same order as those of A and C. On looking for solu- 
tions of the form p = p,e* we obtain 

(1A? -+-e BA— Ag)p = 0, (7) 


where again we have omitted a suffix zero from vector p. We denote the 


elements of the symmetric matrix B by b,;. The condition that equation 


(7) should be consistent now requires that there are 2n possible values of 
A, which are the roots of a polynomial equation with real coefficients. 
Hence, if € is sufficiently small the latent roots occur in complex conjugate 
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pairs, and so the corresponding latent vectors will also occur in conjugate 
pairs. We construct the (n x) matrices A and P from these (exact) roots 
and vectors by selecting one root from each conjugate pair, namely, that 
with a positive imaginary part, together with the corresponding latent 
vector. We then have n latent roots and n vectors with which to form 
the diagonal elements of A and the columns of P respectively. We then 
have PA?+-eBPA—A2P = 0. (8) 

The coefficients of the polynomial equation whose roots are the A; are 
themselves polynomials in e. The A’s may therefore be viewed as analytic 
functions of e, and may be expanded as power series in e. This will also 
be true of the elements of the vectors p;. We may write, therefore, 

Ay = Age teaAyte'Ag,+..., fora = 1, 2,..., n. 

Thus A,,; gives the coefficient of e” in the expansion for A;. Now when 
« = 0 we return to the undamped system, so the element A); gives the 
natural frequency of the (undamped) ith principal mode (Aj; = —w?, as 
defined at the end of section 2). If we define A, to be the diagonal matrix 
whose ith diagonal element is A,,, then we may summarize the above 
equations in the form 


ri? 


A = Ayt+eA,+e?A,+.... (9) 

From equation (8) it is clear that when « = 0 we may choose P to be 
the unit matrix, provided all the roots are distinct. When a root of 
multiplicity « occurs, then a elements of each of the corresponding latent 
vectors are again indeterminate (cf. Appendix I). However, if these 
vectors are linearly independent we can find a transformation which will 
reduce P to the unit matrix without destroying the orthogonality pro- 
perties of the modes. Thus we may take the first approximation to P to 
be the unit matrix and write 

P = 1+eR+eS+.... (10) 
If we also suppose that when e + 0 the columns of P are normalized so 
that the ith element of p; is unity, then the diagonal elements of R&, S.,..., 
are all zero, while the off-diagonal elements are O(1). 

We now substitute the series (9) and (10) into equation (8), and by 
equating the coefficients of successive powers of the parameter € to zero, 
we obtain corresponding approximations to the solutions given by A 
and P. 


The coefficients of « give 
(RA§—Aj R)+2A,A,+ BA, = 0. (11) 


We now consider a diagonal element of this matrix equation (remem- 
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bering that Ay, A, are diagonal, while the diagonal elements of R are all 
zero). We obtain 
1, = —3 


hb;;, since Ay; ~ 0. (12) 
The off-diagonal elements in equation (11) give 


ON fae i tj. (13) 


ij : 4 9 
* (AG; —AG;) 


This equation indicates that our analysis will break down if (Aj; —Aj;) = O(e) 
while the numerator b,; Ay; is non-zero. Thus we must stipulate that 
(AZ; —Ag;) >. for all i, j and i + j, (14) 
if the expansions (9) and (10) are to have any hope of convergence. 
Equations (12) and (13) give the first-order corrections to the solutions 
of equation (7), and agree with results given in (1) and (2). 
The coefficient of «? in the expanded form of equation (8) gives 


(SA2— A2 S)+ (A?+2A, A,)+2RA,A,+ B(A,+ RA,) = 9. (15) 


The (i,7) element of this equation gives 


—_— 


n 
(AE, + 2g; Aes) +844 Are tes S by 71; 0. (16) 
I=1 


Since we are assuming that Ay; 4 0, we use equation (12) and obtain 


bh? n 
u 1 % ard 
Adi 8d 3S bari: (17) 
OA 1 
and the 7; are, of course, given by equation (13). 
The off-diagonal elements of equation (15) give 


nm 


8 5j(Aj — Abi) + 2 j Aoj Ang +555 Ady 4 doy 2 Paty =O (% 


By using equations (12) and (13) we then obtain 


bis Bi (Abi T Aji) Aoj : 


: \ SY bar; (19) 
2(0B— AB)?” RAB) 


Equations (17) and (19) define the matrices A, and S of the expansions 
(9) and (10). These results are still sufficiently simple to be preferable to 
an exact analysis in some cases, particularly if m is large. This will 
probably not be true of the expressions for the next terms in the series 
for A and P. 

Although null latent roots for the unperturbed system have been ex- 
cluded from the above, some of the results obtained can be applied if there 
is one such latent root, A? = 0, of equation (2). In this case a zero appears 
on the diagonal of matrix A,. We suppose, then, that Ay,, = 0, while no 
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other diagonal element of A, is zero. Equation (11) now tells us nothing 


about A,,,, but (13) is still true and gives 


Im: 


b.. 
ng = — x and Tim = 9. (20) 
3 


= —b 
as one might expect, while (19) gives 
bb IX 
mj ~3) 
8inj — am > Bin Vij 
ay i 7-4 
and the 7; are given by equation (13). Finally, by substituting (20) and 


g. 7 jm —_ (23) 
——— r2. . “~e 
0) 
In order to find A,,, we must consider the coefficients of €* in the expan- 
sion of equation (8). If this is done it is found that 


A 


It should be noted that this analysis is applied to only one of a pair of 


Equation (16) then gives Aum =m 


(21) in (18), we get h 


0. 


2m ~ 


null roots of the equation obtained from (7) on putting « == 0. There is 
still one root of equation (7) which is exactly zero, even when e« + 0. 
This becomes apparent if we consider the secular determinant of this 
equation. The root A,, for which estimates have just been found must, in 
fact, be real; so one might guess that A,,,,, A . are all zero. 


m?> “Amr 


4. Some natural frequencies exactly equal 

The natural frequencies referred to are, of course, those of the undamped 
system. The possibility of such a case arising in practice is very remote 
(with one exception to be discussed), but we need to consider it before 
going on to the more likely (and troublesome) case in which some natural 
frequencies are very close together. This problem is of particular interest 
as it belongs to the variety in which a small cause (introduction of damping 
forces) gives rise to a large effect. The physical explanation for this is 
well known. When an oscillation occurs in the ith principal mode, the 
coefficient 8;,p,; in equation (6) represents a slight disturbing force whose 
frequency is very nearly w;. Knowledge of the resonance phenomenon 
then tells us that if w; is very close to w,, then even though £,, is small, 
the resultant motion in the jth mode may ke of large amplitude. 

We will consider the case in which a of the roots of equation (2) are 
exactly equal, while the remaining (n—.«) roots are distinct, and none of 
the roots are zero. We order the latent roots so that the equal roots take 
the first 2 positions. If we now make the physically sensible assumption 
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that there are « linearly independent latent vectors (modes of vibration) 
associated with the multiple root, then it is shown in Appendix I that we 
can still find a matrix Q such that Q’AQ = I and Q’CQ = —A? (as on 
p. 140). But now Q is not unique, for the first « columns of Q (defined as 
matrix Q,) are not defined to within an arbitrary orthogonal transforma- 
tion. Thus, if @ is an orthogonal matrix, and we write Q = [Q,.Q,-.]. 


then the matrix Q* =[0,4, Q,-a] 
will still give a transformation to principal coordinates. On applying such 
a transformation to equation (5) we now obtain 
p+ B*b—Aip = 0, (24) 
G’B,G G'B,, 
Bi,O By | 


and B,,, By, Byz are sub-matrices of B of order («x a), [a x (n—a)], and 


where B* — «B* = eQ*’ BQ* = « 


|(m—a) x (n—«)] respectively. 

Now the difficulty met in section 3 when equal roots occurred is elimi- 
nated if the 6,;; in the numerator of equation (13) are all zero. That is, 
for i,j <a in this case. This desirable state of affairs is attained if we 
choose G to be that orthogonal matrix which reduces @’ B,, G to diagonal 
form. This matrix will be unique if the latent roots of matrix B,, are all 
distinet, which will generally be the case. We denote the elements of B* 
by b*, and with the above choice of G we have 

b*, 
and note that b%, = b,, for i and j > a. 


0, fori,j <aandi +), (26) 


We may now replace equation (8) by 

PA?+«B* PA—A3 P = 0, 27) 
and the expansions (9) and (10) for P and A may still be applied. On 
substituting these into equation (27) and equating coefficients of € to zero 
we obtain (RA2—A2 R)+2A,A,+B* A, = 0. (28) 
Now the first « diagonal elements of A? are equal, so that the leading 
minor of (RA?— A? R) of order (a X x) is null. If we extract the correspond- 
ing leading minors from each term of equation (28) noting the definition 
of equation (25), we get 

2(A,).+G@’ By, G = 0, (29) 

where (A,), is the leading minor of A, of order (a <a). Thus, 

Au = —4bh,, fort <a. (30) 
In practice, the matrix @ and the elements A,; can be constructed from 
the orthonormal latent vectors and the latent roots of the equation 


(B,,+2]p)g = 0. (31) 
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If « = 2, we can evaluate A,, and A,, explicitly as 


Ai Are — }(b,,+-592)+ }[ (6;,;—be2)?+-453,]!. (32) 


From this point we can proceed with the perturbation method as 
described in section 3. It is found that the approximations to the latent 
roots given by equations (12) and (17) are still formally true provided we 
replace b by b* wherever it appears. 

For the elements of matrices R and S we obtain 


by oj 


aah fori > «aorj >a, 


ual 


_ 4 bir, fori <aandj < 


 lea+ 


b¥ ; OF, (AG: + ~Aoy) a 1 : 

: a la LE b* fi . ~ 
| 2(A2; a 2)? ayy > Ail ijs ort x or ) x 

at Oe (34) 
| A rj(2A oA —Ai) 


: - 
No (9; + os) b* —b* ? 2 bi, 8), for i Sa and j < < 


It must also be born in mind that r;; = 8,;; = 0 for i = 1, 2,..., . 


The repeated root which is most likely to occur in practice is in fact 
zer>, for it is quite possible for more than one coordinate of the system 
to be unrestrained. If we suppose that A,; = 0 for i < a, and that the 
roots Ay, are distinct and non-zero for ¢ > a, then it is again found to be 
necessary to reduce the matrix @’ B,, @ of (25) to diagonal form. Having 
done this the results given below can be obtained; the analysis is a little 
more awkward than the above, as it is necessary to go as far as the 
coefficient of «4 in the expansion of (27) in order to complete the matrix S. 
We find that 


Ai = — b*. for i ~ Qa, Ai = ‘ bhi; for i - a, (35) 
Ay, = 0 fori <a, (36) 


and equation (17) (with 6 replaced by 6*) holds for i > 


2 


0 for j <a 
bf; dog 
dg, — 23, 
_ by 
ri, 

=o 


for j > 
for i > a andj < a, 


fori < aandj <a, 
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and equation (19) (with 6 replaced by 6*) holds for 7) > a. Some simpli- 
fication occurs in equations (37) and (19) when i < a and j > a, for then 
Ao = 0. 

It should again be noted that in addition to the roots with non-zero 
corrections given by equation (35), there will still be « roots of the damped 
system which are exactly zero. (Ref. the last paragraph of section 3.) 


5. Numerical example 
We consider the matrices 
A= jf : 4 -—6] and C 
4 —20 
-§6 —20 18 - 18 
The matrices Q and A, as defined before equation (6) are then 
Q . 3 l 1], and A, = 
0 l/v2 2 
l l/jv2 3 
i.e. Q AQ = I and Q’CQ -Af. We notice also that any linear com- 
bination of the first two columns of Q is also a latent vector with latent 
root 2. 

The matrix 8 of equation (5) is chosen to be diagonal. The non-zero 
elements are given by the harmonic means of the corresponding elements 
of matrices A and C' multiplied by e. It is then found that the symmetric 
matrix B of equation (7) is 

0-827 963 0-176 811 0-327 865 
0-741 734 0-067 O80 
5-761 640 


Equation (31) is in this case 


» + 0-827 963) 0-176 811 2" 0. 


0-176 811 (2u+-0-741 734) ] |g. 
The latent roots of this equation give 

Ay, = —0-301 428, Ai. = —0-483 420. 
After normalizing the latent vectors we obtain the orthogonal matrix 

G 0-617 696 0-786 417). 
0-786 417 0-617 peed 

Matrix B* is now given by equation (25): 

B* 0-602 856 0 

0-966 840 
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The new principal modes of vibration for the undamped system are defined 


by 0-061 61 


Q* = 1-223 194 


1 


» 


6 
l 


0-556 68 
—1-173 777 


—0-436 777 2 
—0-349 639 3 


We now apply the results obtained in section 4. The following values 
for the latent roots are obtained when e = 0-01. The last row of the table 
gives the roots calculated exactly to seven decimal places for comparison. 








im Te im 
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©9999 9387 ” 
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000730 165 09999 972 00048 317 ©9999 922 -0'0288 o84 1°4139 























The first of the two columns of complex numbers given below is the 
last column of (/+¢R-+«?S) of equation (10). The second column of this 
table gives the exact latent vector for comparison. 





(1 +h +S) Exact solution 





re im re im 
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6. Some natural frequencies nearly equal 

This case is more difficult than that in which some natural frequencies 
are exactly equal. The reason for this is that the effects of small damping 
forces may still be large as a result of near resonance conditions, but since 
the natural frequencies of the undamped system are, in fact, distinct, the 
matrix ( introduced in section 4 is no longer at our disposal. However, 
the line of attack used in that section appears to be the only one which 
will give a general solution by the perturbation method. We therefore 
proceed to modify equations (2) for the undamped system in order to 
produce a system with exactly equal roots; and then allow two distinct 
perturbations of the system to occur simultaneously. One of these is in- 
tended to correct the natural frequencies, while the other accounts for the 
damping forces as described above. The orthogonal transformation matrix 
@ is then at our disposal once more, and its value is determined by the 
combination of perturbations. 
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We write matrix A appearing in equation (2) as the sum of two sym- 
metric matrices; A = © +A, where 7 is a small scalar number. Equation 
(2) then becomes [((O+-nA)A2+-C]q = 0, (39) 
and we suppose that the first «a roots of this equation are almost equal. 
It is shown in Appendix II that © and A can be evaluated so that the first 
x roots of the above equation are exactly equal when 7» = 0, while the 
remaining (n—.«) roots are identical with those associated with matrices 
A and C. Furthermore, if Q is the matrix of latent vectors corresponding 
to the non-zero value of n, then Q’OQ and Q’AQ are also diagonal matrices; 
the last (n —a) diagonal elements of Q’ AQ (= A) being zero. 

We now normalize the columns of Q so that Q’OQ I. (We assume 
that the varied matrix © obtained from A still corresponds to a positive 
definite quadratic form, q’Qq.) If we now introduce the damping coeffi- 
cients into equation (39) and transform to the principal coordinates given 
by q = Q*p, we obtain as in equation (24), 

(J +7A*)p+e«B*tp—Azp = 0, (40) 


where B* is detined in equation (25), and A* when divided into corre- 


sponding sub-matrices is given by 


| VAG ‘| 


0 0 

We write A* = GA, G, and note that although A, is defined in Appendix 
Il to be a diagonal matrix, A¥ is not necessarily so (i.e. G is not chosen 
in this way.) 

Now the parameters « and 7 are quite independent of each other, but 
in order to proceed further with the problem in which both perturbations 
, occur simultaneously, we restrict our attention to cases in which ¢€ and » 
are small and of the same order of magnitude. If we examine the definition 
of « more closely, we find that in equation (6) we must have 

B;;, <<, for w; #9, 
and Bi; l, ,= 0. (41) 
From the definition of » we require 
w (w?),,, (w?),,. fori < a, (42) 
where (w*),, is the mean value of the squared frequencies corresponding 
to the first a modes. Since we only consider those cases in which (w?),, 
is well removed from zero, our criterion for using this simultaneous varia- 
tion of parameters is that the maximum values taken by 
Bis 7 
(a) =, oF B; ifw,=—0 (i n), (43) 


WwW; 





FREE VIBRATIONS OF LIGHTLY DAMPED SYSTEMS 149 


w? —(w?),,| 


(w*)m 


must be very much less than one, and must both be of the same order 
of smallness. On this assumption we now replace 7A* in equation (40) by 
eV*, and write the elements of V* as é. 

We can write, as in equations (8) and (27): 


and (6) (¢ <a), (44) 


(1-4+-<¥*)PA2+«B*PA—A3 P = 0. (45) 


If we now substitute expansions (9) and (10) in this equation and equate 
the coefficients of successive powers of ¢« to zero, the coefficient of « gives 


2(Ay)a +@(By +V, An )@ = 0, 


and V, is the leading minor of V of order (aa). Thus, we must now 
choose @ to be that orthogonal matrix which will reduce G’( By, +-V,Ag,)@ 
to diagonal form. Although this is in principle no more difficult than the 
reduction met with in section 4, in practice the problem is complicated by 
the fact that the elements of (B,,+V,A9,) are complex numbers. This 
implies that, in general, the elements of @ will also be complex. When G 
is chosen in this way, we obtain 

Aug = —HOR+E% AQ) «fori < a, (46) 
while b¥, +68 Ay, = 0, fori andj <a, andi ¥j. (47) 


We still have, of course, A,; = —4),, for «a <i <n. The complete first 
and second approximations are then given by 


=? . 
Agi = ay, Pee 4 2 OT. (48) 


fori > aorj >a, 


b* oj 

{aaa 
l 

2(yy—Ayy) 


(ryt > bina, fori <aandj <a. 





bE (AR, +AR,)A , ro { es n } 
for i >a or j > a, 
a n 
+a ty 2 Satyt 2 big) 


fort <aandj <a. 


2(A,;—Ary) | Yor 
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7. Numerical example 
We consider the symmetric matrix 
] 33:6596 —4-5208 — 55404 
63-9744 40-5584 —20-2008 
17-9796 
and the matrix C given in the example of section 5. The (exact) latent 
roots A? of equation (2) are, in the notation of Appendix IT, 
£, = —1-02, &, = —0-98, £, = —2-00, 
while the latent vectors are those given by matrix Q of section 5. If we 
choose » = 0-02, and use equation (IV) of Appendix IT we find 
(A? =) X, = [—1-0 0 0 : X, = f—1-0 0 Oi. 
0 —1-0 0 +10 O 
0 0 — 2-0 0 0 0 
We calculate Q’AQ and find the diagonal elements to be (exactly), 
a, = a he oe oes a, = |. 
1-02 _ 0-98 
Then equations (V) give 
6=1, =1, 0,=1; 


. 


i.e. © = J, and it is not necessary to renormalize the columns of Q. In 


fact the matrix © is matrix A of section 5, and this example was con- 
structed so that comparison could be made with that of section 5. 
Equations (VI) give 


We consider the damping forces represented by matrix B of section 5, 
and again consider « = 0-01. The relative magnitudes of the two per- 
turbations as given by (43) and (44) are 

ebg, (0-01) x (5-8) i —(w?),,  1:02—1-00 


: = (0-04 and - ——" = ae sz OOD: 
Wig 1-4 (w*) 





m 
so we may apply the method of section 6. 
Since » = 2€ we obtain 
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The values of A,, and A,, and the columns of G are given by the latent 
roots and vectors of 
2u-+[0-827 963—i(1-960 784)] 0-176 811 9.) =9, 
0-176 811 2u+[0-741 734—i(2-040 mn 
(cf. section 5). We obtain 
Ay, = —0-414 066-+1(0-976 480), Ay. = —0-370 782—i(1-016 496); 
G = [—1-000 980+72(0-000 042) 0-000 960+-7(0-044 295)]. 
| 0-000 960—17(0-044 295) —1-000 980-+4-7(0-000 od 
We now calculate B* from equation (25), and V¥ = G’V, G, whence we 
can apply the results given in equations (46)-(50). 

The tables given below are comparable to those of section 5. The first 
compares the successive approximations to the latent roots with their 
exact values, while the second table compares the last column of 
(J+«R+e&S) with the exact latent vector. 





1 t=2 





re im re im re im 
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(1+eR+eS), Exact vector 
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O'OOO!I 534 —O°0009 257 | O'OOOI 517 — ©0009 245 
| ie) ° 1°0 ° 














The leading entries in the first of these tables may be a little misleading. 
The values of 1-0 for im(Apg;), im(Ap.) are, of course, purposely removed from 
the undamped natural frequencies which are (—é,)' and (—€,)}!, i.e. 
1-0099 510 and 0-9899 499 respectively. 


8. Discussion of results 

In interpreting the results arrived at it must be borne in mind that the 
quantities A), appearing in the analysis are all pure imaginary numbers (or 
zero). For this purpose, therefore, they are better replaced by putting 
Aoi = tw;, where w, is the natural frequency of the ith normal mode of 
vibration when no damping forces are present. We will consider first the 
theory of sections 3 and 4 in which the perturbations correspond to forces 
which are in phase with the velocity of the system. 

Equations (13), (19), (33), and (34) show that the first and second 
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corrections to the normal modes of the system are pure imaginary numbers 
and real numbers, respectively. The corresponding corrections to the latent 
roots given by equations (12) and (17) are successively real and pure 
imaginary. In fact, the form of equations (11), (15), and those obtained 
from coefficients of higher powers of e, indicates that this pattern is re- 
peated if one proceeds to higher approximations. This fact will account 
in part for the high degree of accuracy obtained in the example of section 5. 

The imaginary part of the roots A,; obtained in that section are, in fact, 
correct to third order; while the real part of the latent vector (which is the 
third column of (/+¢R-+-e?S)) in the second table of section 5 is also 
correet to third order. 

When the system is oscillating about a position of equilibrium which is 
critically stable, then at least one of the w,’s is zero. Equations (20) and 
(37) indicate that to first order accuracy we may say that the normal 
modes in which divergence may occur are unaffected by the introduction 
of small damping forces. For each pair of null latent roots of equation (2) 
there will remain one null root for the damped system, while the second 
becomes a small negative number (i.e. a subsidence) given to the second 
order by equation (35). 

We now turn to the problem discussed in section 6 in which « of the 
latent roots are very close together, and conditions (43), (44) are satisfied. 
The artifice used in that case appears to have no physical counterpart. 
The manipulation of the ‘inertia matrix’ might correspond to varying the 
mass of the system at those points (if any) which are nodes of the last 
(n—a) normal modes. Such a variation might bring the first « roots into 
coincidence, but would unfortunately alter the normal modes too. Our 
artifice must therefore be accepted for what it is, without trying to give 
the process some physical significance. 

Since the perturbations introduced in equation (40) are partly in phase 
with the acceleration, and partly with the velocity, the simplifying features 
discussed above are lost. Thus equations (46)-(50) indicate that we get 
complex numbers for each approximation to the new latent roots and 
vectors. The fact that the orthogonal matrix G is complex implies that 
when the first a latent vectors are chosen to fit the perturbation, instead 
of having the new normal modes described as a sum of the type 


> 7:4,c08a,t 
i=] 


which is the case in section 4, we obtain 


> ¥,4; c0s(w, t+-¥,). 
1 
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9. Other applications 

Although the method developed here has been applied to the classical 
problem of light viscous damping, it may obviously be applied to other 
problems, notably to non-linear perturbations for which the exact solutions 
may be much more difficult to obtain. In particular, solutions can very 
easily be obtained if the damping forces are of hysteretic type, i.e. when 
the damping force is still in phase with velocity, but the energy loss per 
cycle in simple harmonic motion is independent of frequency (Bishop (7)). 
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APPENDIX I 


Simultaneous reduction of two matrices 
We consider symmetric matrices A, C of order n which have rank n, and r(< ”) 
respectively. We consider equation (2) of the main text: 
(2A +C)q = 0. (1) 
If a set of n linearly independent latent vectors are found together with the corre- 
sponding latent roots we have, as in equation (3), 
AQA?+CQ = 0, 
and also MVQWVA+Q’C = 0, 
since A and C are symmetric. Premultiply the first of these two equations by Q’) 
and postmultiply the second by Q@. Then 
AAt+C = 0, 
and A244+06 0, 
where 4 Q’AQ and C = Q’CQ. Thus, 
AA® = AtA, (II) 
Now if the latent roots of equation (I) are all distinct this can only be true if A is 
a diagonal matrix. Hence, in this case we may normalize the columns, q;, of matrix 
@Q in such a way that 
QWAQ=1 and Q’CQ = —A*. (III) 
If there is a root of multiplicity a among the latent roots of equation (I) we order 
these roots so that the equal roots take the first a positions. Then since q;, Q»,...,q, 
are each associated with the latent root Aj, we may write 
(ATA + CBG, + 8242+ ---+B.4,) = 9, 
where the 8's are scalar constants; i.e. any linear combination of the first a vectors 
is also a latent vector of equation (I) with latent root Aj. Thus if Q, is the matrix 
whose columns are q;, Qs,---. q,, we can always replace it by Q, /’, where F is a non- 
singular square matrix of order a. Thus, we replace 


Q=(2,,0,.] by QO =[2,F, Qn.]- 
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A Q’AQU — A(@. F. Qu...) i. 2 0 I (IV) 
where A, and A,_, are the obvious partitions of A. Equations (II) and (III) tell 
us that with a suitable normalization 4, _, can be replaced by a unit matrix, while 
1. is not necessarily diagonal. In equation (IV) we now choos F to be one of that 
class of matrices which will transform 4, to the unit matrix. Clearly, if @ is an 
orthogonal matrix of order « and F is a member of the above class, then FG is also 
a member, for if ~I 
F’A _F rs 


(FG) A (FG) G1.G I 


N—a 


Thus, when there are a latent roots of equation (I) equal to each other, we can 
still construct a matrix of latent vectors Q such that 
AQ I and Q’CQ A®; 
but the sub-matrix Q, of Q is not defined to within an orthogonal transformation. 
We note also that the simultaneous reduction of A and C to diagonal form is 


not affected by the presence of latent roots which are zero (i.e. when the rank of C 
is less than nm). 


APPENDIX II 


Construction for a problem with exactly equal latent roots 


We consider the (symmetric) matrix equation 


(AE+C)q = 0, (1) 
where € is a scalar number. We suppose that all the latent roots €;, and latent 
vectors q;, are known; and that the first a of these roots are nearly equal, while the 
remainder are distinct. We are to construct a matrix ©, nearly equal to A, such that 
haa 
the equation (QE C\q 0 (11) 
has a latent roots exactly equal, and the remaining (m—«a) roots in common with 


the last (m—a) roots of equation (I). 


We write the diagonal matrix of latent roots of equation (i) as the sum of two 
diagonal matrices 


X Xot+ 7X1, (III) 
where » is a small scalar number. Or, equating diagonal elements, 


Ei = foi + nbs (IV) 


and we define €,; = 0 fora <i < n; and fori < a we define &; = &,, a mean value 


of &,, &,.... €,. The matrix ¥, constructed in this way is to be the matrix of latent 
roots of (II); i.e. this equation is to have the first a roots equal to &). 

If Q, R are the matrices of latent vectors (not normalized) of (1) and (II), we know 
from Appendix I that we may write 


Q’'AQ A, Q’CR = ¢ - 
and ROR = 0, R’'CR = —X,9, 
where A and © are diagonal matrices. Since Q’CQ and R’CR are both diagonal 


matrices we may suppose that the columns of R are proportional to the corre- 


sponding columns of Q. In fact, since © is still to be determined, we may write 
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R = Q, and use the above relations to determine ©. On this assumption we get 
X,0 = XA, 

which defines ©, and we then have © = (Q’)-'0Q-!, which is clearly symmetric. 

However, in the main text we only require 0, so @ need not be calculated explicitly. 

Let the diagonal elements of 4 and © be a; and 6, respectively (/ 1, 2,...,). We 

now have (using (IV)) 


(Vv) 


{ a; (1 + ne a, fori : 
| i 


i for a - 
Thus, we may write A = O+nA, 
where A is a diagonal matrix with diagonal elements 
8; = — resiitle g (V1) 
0, fora <?t<n. 

We see that the matrix © constructed in this way is nearly equal to A as required, 
but has the further useful property that the latent vectors of (I) and (II) are the 
same. 
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SUMMARY 


It is a common hypothesis that the error involved in a certain numerical process 
is an analytic function of a parameter which specifies in some way the closeness of 
the approximations used. We consider this hypothesis in detail for difference 
approximations to differential equation eigenvalue problems. Our approach is to 
construct the solution of this difference equation by determining a recurrence for 
the terms in its Taylor series. Existence is then guaranteed by the convergence of 
this series. The solution thus obtained is used as the basis for a discussion of methods 
by which h?-extrapolation of the eigenvalues and eigenvectors of the difference 
equation could be obtained. 


1. Introduction+ 

To make our problem definite we restrict ourselves to the eigenvalue 
problems of second-order, self-adjoint differential equations, although our 
method clearly has much wider application. The results could readily be 
duplicated in a particular case by following the procedure outlined in this 
paper. 

We proceed in the usual way, and approximate to the differential 
equation by a difference equation on a finite grid of points of mesh length 
h. As we increase the number of points in the grid and let h + 0, we expect 
the eigenvalues of this difference equation to tend to the eigenvalues of 
the differential equation. We can write this as 

lim A,,(h) =A, 
for the pth eigenvalue. We also ask that the components of the pth eigen- 
vector of the difference equation tend to the values of the pth eigenfunc- 
tion evaluated at the ordinate points x = th. That is 


lim d,,(x = ih,h) = $,(2), 


h-+0 


where here we regard x as a fixed ordinate point and i (an integer) and h 
as variables, and where we set 


$,(@,h) = $,(#) #0 
for all A and fixed % to make definite the sealing of the eigenvectors. 


+ The argument in sections 1-3 follows closely that of Fox (1), 332-42. 


{[Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 2, 1960) 
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This behaviour is certainly guaranteed if we can form the expansions 


a hm d™x 0 
Mh) = > = 


<— m! ~dh™ 
m=0 


x 


_ al = pm/ - an. me! OMP(z = ih, 0)| 
and = ¢(z = jh, h) = - — > G-i a ae 
m= .rTr+e m 


In all subsequent discussion we will write 


empl = ih,O) — d* {a"dx = ih,O)\ dad (x) 


oh" J dx* @h 


 aatohe dx 
to emphasize that the coefficients of h” in (2) are functions of x alone. 
Note that we could just as well have written (2) as 


‘ “ hm : 
d(x = th,h) = bs — bm(t = th). 
P 


m 


However, convergence of this expansion gives to ¢,, the meaning €"¢/éh"; 
and our notation anticipates this result. 
The most general differential equation of the type we consider can be 


written, after suitable transformation, in the form 


{2 
> r(x)\A—q(x)}p = 0. (3) 


Together with (3) we consider the class of boundary conditions 


Ay $(0)4 Ree =0 


and A, 4(1)+B, me 0 (5) 


given at the ends of the interval (0,1). We put r(z) = 1 in our subsequent 
calculations as this simplifies the arithmetic—it is not difficult to return 
to the more general case. 


2. Approximation of the boundary conditions 


We first consider what restrictions central difference approximations to 
(4) and (5) place on the expansion (2). The simplest representation of 


(4) is 
) is By, 


Ay d(0, h) + 2h t 


d(h,h)—d(—h,h)} = 0. (6) 


This must be an identity in h. Therefore, when (2) is substituted in (6), 
the coefficients of powers of A must vanish. This leads to a series of 
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relations, the first few of which are 
h®: A, d(0)+-B eH) _ 


0 dx 
ai: 4,209), p & 29(0) _ 


och © dx ch 


pes 45280), yp {¢ 280) , 1240) " 
ie %\dx ch? 3. dz? | : 


2, 


It is to be noted that the first two relations are similar to the given 
boundary condition (4); but that this is no longer true for the third. If 
we take a more accurate representation of (4)—ignoring fifth differences— 
we get 

ae 


A,.¢(0,h) 4 / 
0 PN) Th 


(2h, h)+-8d(h, h) —8d( —h, hh) +46(—2h,h)} = 0. 


(10) 
This fails to reproduce the boundary condition in the coefficient of h* 
giving 


ph: 4 OO phe OO SEEN (11) 
ch! "\dx ch! 5 dz | 
In general our ability to reproduce the given boundary condition in the 
coefficient of h” depends on the accuracy of our difference approximations. 
It would seem from the above examples that central difference formulae 
including differences up to the (2m—1)th fail in the coefficient of h?”. 


3. Formal recurrence for the solution of the difference equation 


The simplest difference equation approximating to (3) is obtained by 


setting “ d2 i 
dx? 
and ignoring the higher order differences. This gives the equation 
d{ (i +- lh, h) —2d(ih, h) +-d((i— Ih, hh) +h2(A(A)—q(th) )d(th,h) = 0, (12) 
which can be written as the matrix eigenvalue problem 
[A+h?A(h) Tip = 0. 


If the boundary condition (6) has been used, then the matrix A is 


~2—h*q(0)+2hA,/ B, 2 
| 2—h7g(h) 1... 


| : | 


The eigenvector @ has the components ¢; = ¢(th, h). 
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Substituting (1) and (2) into (12) and writing th = x to indicate an 
arbitrary ordinate point we get, on equating to zero the coefficients of 


powers of h, the series of equations the first few members of which are 
given below: 


[a2 


\ da? 


h?: +-(A -qix))\¢ — (), (13) 


3. (ala OF — _2 
ai Maia = — qn ee) (14) 


eb _ d?r 9 fdr cd(x) 1 d*p(x) 


ha: 4 ~(x : . ry—e = i 15 
: 2 I | ch? ane?) dh ch 6 dx* as 


\%h dA 4(2) _g PA Ed(x) 
\ch? —s dh3" "dh? Gh 


q(x)) 


dr h(x) I d* &d(x) 
"dh oh? 2dzx* ch’ 
\etd d4r dX &d(x) 


UN ape = — ae —4 is oh 


6 d*r ed(r) Pa Ad(x) d* d(x) 1 d®p(x) 


‘dh? Gh? dh @h® ~~ dxt ch? 15 dx® 
The general equation is 
, ]? je"d : 
mee: a | _ : ass . 
t (a (A 9)) Zim V(x) 


m—1 


' m—p) ap 
where V(2) = — 7 m! d™-P) Erd 


7 p\(m—p)! dh™-» eh» 


p=0 
[4(m + 2)] 4 P 
d2r am+2- 2rd m' 


da Ehm2-® Fim 42 —I)!" 

We see that (13) is the original differential equation (3), and that in all 
subsequent equations we regain the left-hand side of (3) operating on 
emd(a) ch™:; but now forcing terms appear on the right. 


4. Solution of this recurrence 


We now discuss the way in which the terms in expansions (1) and (2) 
are determined by the manner in which they appear in the above equations. 
In equation (14) we know that é¢(x)/ch always satisfies the same boundary 
conditions as the eigenfunction ¢(x) of the original equation. Therefore 
the corresponding A is again an eigenvalue; and, for a solution to be 
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possible, the right-hand side must be orthogonal to ¢(x) over the range 
(0,1) of x. That is 1 


dd 
_d? dr — 0 
[ai salen, 


whence ag 0 
dh 
We thus regain the original equation, and hence ¢¢(x)/ch is ¢(x) apart 
from a numerical factor. This factor is zero for we agreed to fix ¢(2 = Z, h) 
as a constant (non-zero and independent of h) for all h, and hence 
o"'d(x z)/ch™ 0 (m > 1). 
The general equation (18) not only has the right-hand side non-zero, 


. 


0 


but also satisfies the inhomogeneous boundary conditions 


' emd(0) B d emd(0) (0), 


ov eh™ "dx Gh™ sa 


F omd(1) RB d omd(1) 


and y. 
1 chm 1dx @h™ 


> H,,(1). 


The condition that (18) has a solution can be expressed in the form (see 
(2) 214-16) 


m 


G(O)dd(0)/dx}  H, (1){dp(1)/dar} 


(19) 


1 
[ Af tte ; x} d(x) dx 
0 


Applying (19) to equation (15) we get 

B, d34(0) 
3 dx 
B, d*4(1) 
3 a 


G,(0) 


H,(0) 


and 
1 1 
d*\ ” 4 Py B, d4(0) dd(0) BB, d¥g(1) dd(1) 1 f d*d 
dh? . “34, dx? dx 3A, dx® dx 6 | dx* 


0 0 


d dx. 


This simplifies to 


1 1 
dh {1 d34(1) 1 d34(0) Ifa, 2) f os 
=m . : a 
dh? \3 dx * 3 dx #(0) 6. aa? ™|/ | — 
0 0 


Equation (19) can always be solved for d”A/dh™ because we can assume 
that, at this stage, all other terms and their derivatives with respect to 
x are known quantities. The solution of (18) is determinate to within 
an arbitrary multiple of the original eigenfunction ¢. This is fixed by the 
condition é"d(z)/ch™ 0. 
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We now show that "¢/ch™ and d™A/dh™ are uniquely determined by 
the condition that (18) should have a solution satisfying the imposed 
boundary conditions together with the condition ¢"¢(Z)/ch™ = 0. 

We prove this by induction. It certainly holds for m = 2 for, if we 
denote the difference between two solutions satisfying the above condi- 
tions by A(@*d/ch?), we have the ety differential equation satisfied 

(@* . y_@)\ -(45, 
: q) a? : A 
laze 9) ~e 
also the ori,sinal boundary conditions (4) and (5), and the adjoined condi- 
tion A(é*A(z) ch?) = 0, whence 


d*X e*p(x) 
-=0, and A $(: -= 0. 
dh? ch? 
Now if we assume that the result is true for c'¢ ch’ (¢ < m), we can repeat 
the above argument exactly for e"d/ch™ and d”™A dh”, and hence the 
uniqueness follows. 


5. Convergence of the formal series 
We can now establish the convergence of expansions (1) and (2) under 
the assumption that g(x) is a regular function with a non-zero radius of 
convergence p about every point x in (0,1). We know from standard 
theorems that ¢ and cd ch” are regular functions with the same radius 
of convergence as g(x). Let the minimum value of p for x in (0,1) be not 
less than p. 
We have the estimate 
| dP og| “3 p! K. 
\dax? che ~ pP 
where A is independent of p, g, and x and is bounded in (0, 1) since 
od ch® has no singularities in this range, and vanishes at 4 = 7. 
We can now form an expansion which dominates (2). We have that the 
modulus of the mth term in (2) is less than 


a — 4 m -_ 
' fam : ts Jie J =) +) a. 
al aa" (m—p)! | p—t™-? p™ SF ad 


h™ m'K om l =m—p Rete ee «A l =p 
jimh™m S RE his h K > | - 


-t 


c= |\j—t mi ~ Kexp —_ P_\ 
: 5 \ 


The series 


converges for |j—il\h <p p- 
We proceed to prove similarly the convergence of (1). To estimate 


d™\/dh™\| we evaluate (18) at x = . We know that the left-hand side 
5092 .50 M 
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is bounded—in fact €"d(Z)/ch™ = 0, and we can take (d?/dx?)(@"4¢(Z)/éh™) 
into the second sum on the right. Writing K, = K/¢(#) we have, since 
all terms in the first sum vanish except (d’A/dh™)d(2), 


[i(m + 2)} , [i(m+2)] 49-5 
d™r S (2p)! ,. m! 2K,m! on pm +3—Sp) 


dh™ — p* '(2p)! (m4 2—2p)! | am — (m4 2- 2p)! 


p=1 
2K,m!{ sinhp: m odd 
p”*? | coshp: m even, 
2K,m! - 
; ~- cosh p. 


and hence —— 
am 


, a p* p™ . 


m=O 


So 2K, cosh p hm 


The series 


which dominates the series (1) for A(A), converges for h < p. We have 
thus proved that expansions (1) and (2) exist, and that they satisfy (12) 
for h - Pe 

6. Proof of the ‘deferred approach to the limit’ 

The result that forms the basis of Richardson’s technique of deferred 
approach to the limit can now be deduced. 

If central difference formulae are used to approximate to the boundary 
conditions (4) and (5) of the differential equation (3), then all odd deriva- 
tives of ¢ and A with respect to h are identically zero. To prove this we 
need the following results. 

(1) The inhomogeneous terms in the boundary conditions on 

{gam + 14(x)} (eh2m+1) 


involve odd derivatives with respect to x of odd derivatives with respect to h 
of degree less than 2m-+-1. 

This follows because (a) the 1/4 term multiplying the difference operator 
means that derivatives of total order 2m+2 are involved; and (6) the 
coefficient of d(2,h) is the negative of the coefficient of 4(—z,h) in the 
central difference expansion of dé/dx about 2 = 0, and a similar result 
holds about « = 1. Thus the even derivatives with respect to x cancel. 


(2) The forcing term on the right-hand side of the differential equation for 
e2™+14h(x)}/(Ch2™*!) consists of odd derivatives of A (db) times even derivatives of 
d (A)—the derivatives being taken with respect to h, and even derivatives 
with respect to x of odd derivatives of & with respect to h. 

With this information we can form an induction. 

If we assume that odd derivatives of ¢ and A with respect to h of 
order up to 2m—1 vanish, then é?"+!¢/eh?"+! must satisfy the boundary 
conditions exactly (@,,,., = H,,,., = 0); and the right-hand side of the 


m 
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differential equation must be equal to —(d?"+1A/dh?™*')¢. We have 
already proved that éd(x)/@h = dA/dh = 0; and, as the circumstances here 
are exactly similar, our result is established. 

This vanishing of the odd derivatives of ¢ and A with respect to h does 
not occur in general if forward or backward differences are used to approxi- 
mate to the boundary conditions. For instance, if we approximate to the 
condition dd¢(0)/dz = 0 by the difference relation Ad(0,h) = 0, the 
boundary condition on @¢/¢h is 


d ad(0) 1 d*d(0). 
dx ch — 2 dz’ 
and we have G, 4 0. This leads to a non-trivial relation for dA/dh; and 
hence we get h-extrapolation only. 


7. h?-Extrapolation of the eigenvalues 
Our equation (20) can be used as the basis for a method of obtaining 
h?-extrapolation of the eigenvalues. We have 


1 


1 
dr  { 1 d34(0) 1 d34(1) _ 1 f dtd | [ . 
dx H0)+ a1) § finde) tied 


dh? | 3 3 dx 6, 
0 


oe  . do, 
za ~-O-8 tee 


d*h > »,,,db ” 
dx' a (A —q) ¢ <q dx +q ¢. 


1 
and | leq 4 "4 dx = q'¢? .° 
0 

so that (20) involves only quantities which lie in the range (0,1), and 
which are known to within O(h?) when the simple difference equation (12) 
has been solved. Thus the error in using the results of (12) and, say, the 
trapezoidal rule to evaluate h?/2! x (20) is O(h*); and we have an h?-extra- 
polation. . 

There is a very close relationship between formula (20) and the difference 
correction given by Fox (1). We show this for the case when boundary 
conditions of the form (6) have been used; and we restrict our attention to 
the integral terms of (20). 

The difference correction AA, with boundary terms suitably interpreted, 
to be applied to the eigenvalues of (12) is just —A?/2! x (20) evaluated by 
the trapezoidal rule with error O(h*). 
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We use the fact that 
ha did(x ih) 
: dat = 
Applying the trapezoidal rule to (20) we get 


h2d2r A h—2 [<< 
4 . >* siny!’ Stg(ih)| / > r {h(ih)}?+ boundary terms +-O(h*). 
y dh? 4 | 1? | / font . 


0 


a. 
= 54d¢(ih)— hh ® a +O(h®). 


m 


We use the notation hg to indicate that we take } the first and last 
members of the sum. 


The difference correction is given by 


PO lina. 

Pp"? 

where m and @* are the right and left eigenvectors of (12) respectively, 

and C is the correction matrix. Here C is the matrix representing the 

operation (h-?/12)5* on @. We see from the matrix form of (12) that, if 

d* = dbo, then d¥ = d; (i = 1,..., m—1), and ¢* = }¢4,,. This establishes 

the required result. The role of the boundary values in difference correc- 
tion is discussed in section 9. 


+- boundary terms, 


8. /°-Extrapolation of the eigenvectors 
We can obtain an /?-extrapolation of the eigenvector by solving (15) 
numerically. In doing this we must use ¢74(#)/eh? = 0 as one of our 
boundary conditions to avoid trouble from the eigensolution of the homo- 
geneous equation. For example, if we assume that difference conditions 
of the form (6) have been used, and that we may fix the scale by setting 
F = 0, d(0,h) = 1 (say), then we get a marching problem for ¢7¢(x)/ch? 
with the initial conditions 
e*4(0) 
oh? 
d ( a) 1 d34(0) 


dx\ ch? 3 dz 


(a) 0, and 


(b) 


and the terminal condition (which gives a numerical check) 
PU oe {d ed(1) | d4(1)) : 
ch? ‘\dx ch? '3 da? | 

In practice one would probably need to repeat the computations at 
least once using the improved value at each stage to reduce the residual 
error. If the marching problem is unstable we could solve (15) as a two- 
point boundary problem using one or other of the initial conditions as 
a check. If the original boundary condition is 4(0) = 0, then we must 
have & ~ 0, and we get two two-point boundary value problems. The 
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condition that (d/dx){e%d(x)/eh?} be continuous at 2 = £ supplies the 
check. 
We could now proceed to obtain an h4-improvement using equation (17) 
provided we have worked to O(h*) in solving for d?A/dh? and 6*¢(x)/éh?. 
As a numerical example consider the equation 
12 
: an (A—a)¢ = 0 (21) 
dx? 
with the boundary conditions 
$(0) + 


_dd(0) = 
dx 


and (1) 4270) = 0. 
dx 


We give below the result of the numerical computation of the first 
eigenvalue; formula (20) was used to compute d?A/dh?. 
g I 
h | A(h) d?X/dh? improved A 
O-1 | —0-82621 0-040 ~0-82601 





The correct value to six places of decimals is —0-825980, 

As a check on the numerical solution of (15) we used inverse interpola- 
tion to determine d?A/dh? such that the terminal boundary condition was 
satisfied. This gave d?A/dh? = —0-041 for h = 0-2. 


9. A method of difference correction 

The process outlined in section 8 suffers from the important disad- 
vantage that the correction terms are defined by differential equations. 
The numerical solution of these (and the construction of the terms appear- 
ing in them) provides another source of error; and a revision procedure 
involving all terms previously calculated is necessary when proceeding to 
a higher order of accuracy. Our work does, however, suggest a modifica- 
tion of the difference correction method described in (1), and we give a 
sketch of this below. 

The exact eigenfunction ¢(2) satisfies the difference equation 

h-*8*h, +- (A—q,)b; = h-*{ 3,84 — 98° +- ...}d; 

at the ordinate point th. 

We write the approximation given by equation (12) as 

h-*527g, + (A—q,)¢; = 9; 
and db; = $;+hC,, 
A= A+h2p, 


using the results of section 6. 
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Now h?('; and h®y are given by the difference equation 
hh} + (A—q, (AC; h2 ud; —h®ph?C;+ 
az h-*{ 304 WOO +..}D; ra 7130" . wo” t --.}h% 
We define h?C'; by the equation 
h- 2h, + (A—q, {hC;! h2jid, +h-*{ 1,84 33% + ...3,. 
As the ignored terms are O(h*), we have that 
WC, = WC, +-h8D,, 
and hu = h®f+hy, 
where /?i is detined by the usual orthogonality condition (that is the 
algebraic condition analogous to equation (19)). 
Now A‘), h'v satisfy the equation 
h *8*{hA D+ (A—q, {kA D} = —hivd, hhc, -h*ih* D, havh?C, — 
h*vh*D,+h-*{ 184— 288+... 2, +h 184 — 186 +... Mat D.. 
Proceeding as before we define h4D, and h4i by the equation 
h-*8*{h8D,\ + (A—q,){h8D,| = —htnd, —h2ah®C, +h- 84 — 388+... Mec, 
ignoring terms which are O(h®), 
This process can obviously be continued; and we have finally 
d; = &+h2C,+h8D,+..., 
A = A+-h2i+h4o+.... (23) 
It would be satistving if at this stage we could make the identification 
1 d*A(h) 
2! dh? © 


and so on; but these relations do not hold. We can, however, say that 


1 ediihh 
( oe, al j= 


2! ch? 


(22) and (23) are asymptotic expansions in the sense that, if we truncate 
these expansions at their nth terms, the remainders will tend to zero like 
h?"; for it can be shown that 


R (d) : kK, S 


L \e2d(ar, A)\ pap 


<~ 2p! oh®? 
pon 


4 l aera) pep. 
: 2p!| dh?P | 


and R,(A) K 
where AK, and Kk, are bounded for a fixed n. 
Also, if the scale has been fixed by setting ¢, = 4,, xh = #, then we 
can adjoin the conditions 
( - D., 


67" h( ah, h) 07"d(ah, 0) ( 


chen F yen 
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Thus the procedure outlined in section 8 can be paralleled for the 
sequence of difference equations. The boundary conditions on the differ- 
ence equations can be derived in the same way as we derived the equations 
themselves. 

For example, we have that 


Ay dby + Byh-pdd, = Byh- {hyd 3 —hyd*...}bo, 
and Ao dot Boh 'nddb_ = 0. 
Hence, 
A g{h®Cy' + By hd h2Cy} = By h- hyd? — syd... }b9+ 
+ Buh 4yd3— hyd? +... 2C,; 

and Af h?Cy} + Boh pd{h2C,} = By h-' hyd 3— dyd'+...}ho.... 

All significant differences can be calculated provided we can extra- 
polate d,; far enough (see (1)). 


10. Another application of our results 

As a final example let us ask when the approximate difference equation 
is satisfied exactly by point values of the eigenfunctions of the dif- 
ferential equation. This means that ¢?"d¢(x)/ch?" — 0 for m > 1, which 
implies that (18) has a solution proportional to 4, so that the condition 
c2"h(&)/eh?" = O ensures its vanishing. This means that the orthogonality 
condition (19) must be trivially satisfied, and hence that d?”"A/dh?” can 
be chosen to make JV,,,(”) vanish identically. Thus we get the differential 
equation 2 
We have G,,, = H,,, = 0, which give as the possible boundary conditions 
d = 0, and dd/dx = 0. 

In these cases it is possible to express A(h) in a closed form. We have 


2m » Ip! 
amy _— om: ( = yn +1ym+i 


dh?” (2m+-2)!' 


A(h) = 2 > ‘ — msi frm 


2m-+-2)! 


m=0 
= a() —cos(vAh)) 
“ th-*sine(S) 


This result is readily verified. Equation (12) becomes 
$541 —(2—hA) gb, +9; = 9. 
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This difference equation has the particular solutions ¢, = e'7®, e-‘8, where 
B is defined by 
(2—h*A) = 2cosB, or A= 4h-*sin? }P. 

From the boundary conditions we get the relations 

(a) (0) = (1) = 0: sinmB 

(b) (0) '(1) 0: sinmB 

(c) (0) = 4'(1) = 0: cosmp 
where | ‘<< m, and h lim. 
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SUMMARY 


An investigation is made of the steady compressible magneto-hydrodynamic 
equations of motion in which certain idealizations are made, namely, that the 
kinematic viscosity is zero, the thermal conductivity is zero, and the magnetic 
diffusivity is zero, These assumptions imply that there is no diffusion of vorticity, 
heat, or magnetic intensity. Under these circumstances it is shown that a Bernoulli 
equation exists containing terms which depend upon the magnetic field in addition 
to the well-known terms which are present in the non-magnetic case. An extended 
vorticity equation is also derived which shows that there exist in general certain 
stream surfaces within a gas at which the magnetic field and the velocity field have 
infinities. As these singular surfaces are approached, the velocity and magnetic 
intensity vectors tend to become parallel. A particular example of a helical flow 
field is discussed which shows that there may be non-uniqueness in certain cases if 
the cylindrical boundaries are not chosen appropriately. 


1. Introduction 

Ix this paper it is proposed to investigate the idealized magneto-hydro- 
dynamic equations for an ionized, compressible gas in which it is assumed 
that molecular viscosity and thermometric conductivity are zero while 
the electrical conductivity is infinite. The vanishing of the kinematic 
viscosity for a non-magnetic fluid implies that the vorticity is conserved 
following the motion of a fluid element; the presence of such viscosity 


will cause the vorticity to diffuse into the gas. Similarly for a non- 


magnetic fluid the vanishing of the thermometric conductivity implies 
that there is no heat diffusion taking place or that an element of gas 
moves with constant entropy. Finally, if the electrical conductivity is 
infinite (or magnetic diffusivity zero) then there will be no diffusion of 
the current filaments into the gas. Accordingly, when solutions of the 
magneto-hydrodynamiec equations are obtained subject to the three condi- 
tions imposed above one can, as far as free motion problems are concerned, 
predict the true behaviour by considering the given fields of vorticity, 
heat, and current to diffuse down their gradients at appropriate rates. 


[Quart. Journ. Mech. and Applied Math., Vol. X1i!, Pt. 2, 1960] 
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We here consider exact solutions of the steady axially symmetric mag- 
neto-hydrodynamic equations and show that the complete three-dimen- 
sional velocity (V), magnetic induction (B), and current (j) fields are 
expressible in terms of two functions: ys, the Stokes stream function in 
compressible flow, and p, the density. These two functions % and p are 
shown to be connected by two equations. The first, a Bernoulli equation 
of an extended kind, contains contributions from the magnetic field in 
addition to the usual terms, the second, a vorticity equation also contains 
contributions from the magnetic field. One of the most interesting features 
of this solution is the presence of 2 singularity in the velocity and magnetic 
induction fields which occurs when 


py’? (yb) 4rr, 


(ys) being ‘the magnetic induction stream function’ and is such that when 
the equality is satisfied, V and B will be infinitely large on this surface 
in three-dimensional space. As this singular surface is approached the 
vectors V and B tend to become parallel. 

\n exact solution is investigated which illustrates the effect of this 
singularity and it is found that a steady solution is unique only in a certain 
restricted region of space, in other circumstances the solution for the 
density is not unique and this may imply that magneto-hydrodynamic 
shocks occur in the corresponding physical situation. Many more special 
cases will have to be investigated before a complete understanding of the 
solution is achieved. 


2. The reduction of the magneto-hydrodynamic equations 

Solutions of the non-viscous hydrodynamic equations of motion when 
the flow is steady and symmetric about an axis are well known; here it 
is proposed to show how this solution can be extended to the hydro- 
magnetic equations. Accordingly, using the usual cylindrical coordinates 
(r, @,z), it will be assumed that the flow and electromagnetic field variables 
are independent of the angle @ and of the time, in which case the complete 
set of equations will be as follows: 


cr 


pus 4 woe 2 Ps” (in B.—j_ Bo) 


oui" = (j. B, 
: 


cr 


Bis 
Ss ALL Bo 


C2 
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be (pru) + = (prw) = 0, 
er C2 


U ¢ Ps 4w~ P = (0, 
or\pY 02z\ pY 


~, = curl B, 
- 


lé ob, 
div B . — (rB,) 4+— = 0, 
ror cz 
E+(V»-B) = 0. curlE = 0. (1.8) 
These equations are valid for a gas whose electrical conductivity is infinite, 
whose thermal conductivity and molecular viscosity are zero, and in which 
the velocity is small compared with the speed of light. 
It follows immediately from (1.4) that we may introduce the Stokes 
stream function &% such that 
an fay 
pru af prw — (1.9) 


and hence from (1.5) it follows that 


p = O(yb)p”, (1.10) 
where @() is a function of % only and may be called the entropy or 
potential temperature function. From (1.8) it is possible to deduce that 
V~B Vé where ¢ is a scalar function and we make the assumption 
that ¢ is a function of r and z only. This makes EF, identically zero; if 
EK, is non-vanishing it implies that the fields are unsteady (1). Accordingly 
we have 

vB,—wBy mJ (1.11) 
cr 
wB,—ub, , (1.12) 
sf,—vB, < ——. (1.13) 
From (1.9) and (1.12) it then follows that 
hep yp _ Aap 


rp cz " 


B 


r 


(1.14) 
rp or 

where A is some function of r and z. If we substitute for u, w, B,, and B, 
in equations (1.11) and (1.13) it follows that 


L ow— By) ® ba : (Av mate. — (1.15) 
pr c C2 C2 


Y or pr 


~ 


and hence d = K(p), (1.16) 


+ Equation (1.5) becomes p = constant = p, in the incompressible case. 
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where A(s) is an arbitrary function of %. It is clear from the above that 
K (is) is the potential function for the electric intensity field E and we have 


E Kw). ve + Sal, (1.17) 


so that E is directed along the normal to the surface % = constant and 
is accordingly perpendicular to both the vectors V and B which clearly 
lie in the tangent plane to the surface ~ = constant which is a surface of 
revolution. 
It follows from (1.15) and (1.16) that 
Av— Bo prk' (ys). 
If we substitute B, and B, from (1.14) in (1.7) it follows that 


dpb) 
Pv de’ 


where (ys) is an arbitrary function of y; it is clear that p(y) is acting as 


pee’ (ip) (1.19) 


i ‘stream function’ for the meridional magnetic induction field, since (1.14) 
can be rewritten in the form 
rB, Culp) rB, Culp) 


Cz cr 


(1.20) 


although it is as readily interpreted as a vector potential component. 

Equation (1.18) is one relation between v and Bg, a second relation can 

be obtained from the second equation of motion (1.2) as follows. 
Equation (1.6) for the current vector j becomes, in the present case, 


wel c Bo. ( B, 6 l 


“ (r By), (1.21) 


c Cz cr ro 


and it then easily follows that the equation of motion (1.2) can be written 
in the form 
; 
plu “ (er) tw (vr) | 1B od (r By) + B, : (rB,)\. 
\ Or C2 | 477| cr " OZ | 

Accordingly when we substitute for u, w, B,, and B, in terms of % we 
obtain 

ow ¢ (vr) 4 ou “(ur een sj > C -(r r Bp) + x = 


Cz Cr Cr C2 rl \ 


Cc 
“4 By) | 
where we have also made use of (1.19); since this equation can be written 
in the form 
OF 2 (ur) 9 2 (op) — 2h), p | _& Z(H), 3), 


er €2z 2 Or cr €z\ 42 ‘| éz Or\ 4a "i 
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it follows that we can obtain its first integral in the form 


ort), B — vy), (1.22) 
4n 


where ¥ (ys) is a function of % only. This is a second relation between v 
and By, and, combined with (1.18) which we can write in the form 


pp (b)u— By = —prk'(p), (1.23) 


we can now solve explicitly for v and is The solutions for v and Bg are 
vr = V(b) + pre WK) (1.2 

and Barf = pp'(b)I be +-pr®K' (yp), (1.2% 

where Bf = *(w). (1.26) 


It is permissible to solve (1.22) and (1.23) for v and By only when 9 is 
not zero and this we shall assume at present. 
We consider now the equations of motion (1.1) and (1.3). In the redue- 
tion of these equations it is convenient to write 
cu ow 
7 pr n’, 
oz oor 
cB, cB, 
—!—-—~ = yf; (1.28) 
cz cr 
7 is the meridional component of the vorticity and y is clearly related to 
the zonal current jg. If we write q for the resultant velocity we then 
obtain from (1.1) and (1.3) the following 


“(r rB,)\, 


Pp (49?) ania ov _ Ps l lry = Fa , 


or ror cr 4n\'* 
ol-< (4q2) —ruy a ans, wee. 2 
\az °° éz{ 


C2 


< ag? )+ Vie" '(b)p”- ee % 2 XP a| 
r\ “p 4np| 


ia 


= alr —E <(B a 


é 2 alle 40 xe 
“) +z Zz — (49 ae <i (yb)py- +1 | 


- = 5° r= (v or) ee for By)} 
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that is, 


o{,. yp | ! , ia xe’ | 
hq? + 7 ©'(yp)py- + 
ér\*4 (y—1)p) y—1 Per" + 4p} 


rBo © B,)\, (1.29) 
: 


is 
a ual Amp ¢ \’ 


rele 
, YP \ HF yyy a 0 xP | 
(y 1 )p) | y—1 p pp 4rp| 


rBg c 


Mh oe ger) “ (rBg)|. (1.30) 
Cz dnp Cz J 


r2| 
The derivation of the Bernoulli equation can be effected most simply as 
follows.t We multiply (1.29) by uw, (1.30) by w and add, then it follows that 


C a lie yp | v a 7 ”s =t.( ae m. 3 
"5. = alt (y—1)p} (us ? 3} 4zrp " or . ez (r Bo). 
(1.31) 


Using (1.22) we can rewrite the right-hand side of this equation in the 


"(u er w : \( wr) Bg (x +wW . Je Bo), 
r\ or rd 4n 4nrp\ or 44 


Ce 


“ 


form 


since the terms arising from ¥ (x) cancel. This expression may be written 


je'(b)v Be | 


( Q) 
e—--w r Bo), 
| 4ezr in| or oa)! @) 


and it then follows from (1.23) that (1.31) can be expressed in the form 


a mus pw \rBe 
(y—1)p} 4n cr C2 


(u : ~ . ) YP + hq? + k ()B,| = (). 

or Cz \(y—l)p 4n j 

Accordingly we can express Bernoulli's equation in the form 
5 K’ ‘ 
YF _ 4 Aq? 4 (Y) -B, F(), (1.32) 

(y—l1)p 4n 
where F (i) is a function of ~ only. 
To determine the vorticity equation we substitute for p from (1.32) into 
(1.29), in which case we obtain 
K'(p) | eby Ob) yn xe’) 


Oln 
fk *Bo' 4 
orl (p) ss r Oy ér | y | Pp p 4rp| 


: , 2 
al" 5, 4p an? 0) 


+ lL am indebted to a referee for this derivation which is simpler and more elegant than 
my own. 
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eff OM) 0 XBL py) AW 
er | y 1?” p ro a (¥) 4n ‘ 


= (Y) € “ (r By) 4 - (vr) ~ (rBg). (1.33) 
4n cor 4nrp cr 


If we substitute for A’(%) from (1.23) on the right-hand side of (1.33) it 


becomes oe 
Xf © (yp) HW) § « (Bo), 
rler 4n ’ 


which we may write in the form 


. rBo , 
hy) +e utwy |. 
r| 4a jor 
Accordingly (1.33) can be written in the form G(é%/er) = 0 and similarly 
(1.30) can be written in the form G(és/cz) = 0; hence 

ff 


G = OM gs 1H py) Oral + 
r 


ye +E HW] =" 


(1.34) 
The equation G = 0 is the vorticity equation for the axially symmetric 
flow. Equations (1.32) and (1.34) are two equations between y% and p and 
there are clearly various forms in which they may be written by using 
(1.24) and (1.25). In general, however, the elimination of either of these 
variables in order to establish one equation in one unknown is not possible 
except under very special circumstances. Since 


rBy = pier” “()} 


a 


the Bernoulli equation can be written in the form 


; = K'() , —¥ - -F 
i= apt + w'() | yur (h)} = Fp): 


hence the function 
K'(p) . ¥ (p)K'(p) 
vr = F : : 
tp TH+ iy) w+ — 0 
also remains constant along a stream surface. 


The vorticity equation (1.34) can be rewritten in the following form, 
which shows that the singularity at o = 0 plays a role of considerable 


(1.35) 
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importance in this problem. We have 


xeon) nw’ (cB, eB, Ljcu  cw| 

imp p = 4npr| cz er| prlez er 
a it & Cus of a} pe’ if =) a(n I) 
pr \ez el er\pror}|  ampr\éz\r ez] | er\r er}} 


l | « (=>) : c *\\ : Me” oy? ' 2): 
pr \e z\pr cz cr\ pr or | 4zpr- 
hence (1.34) can be written in the form 


1 {¢ (==) _! eh 


pr\cz\pr er} — er\pr er} } 


Ob), K" 


- ee ‘Bo » a 
f (us) r p ro rBg “| (ds) | Pad () inp ws). 


é 
(1.36) 
It will be observed that the highest differentials in J have as coefiicient 
the term .o7 p*r*; hence .o = 0 isa singularity for the differential equation 
(1.36). The singularity at .o7 — 0 is clearly so fundamental to the idealized 
magneto-hydrodynamic equations that some relatively simple physical 
expression of the singularity should be possible. This may be shown as 
follows. The equations (1.9) and (1.20) lead to 
l 
py’ (ab) 
v ¥ (ab) +(1 4r)pr2p' (bh) AK’ (pb) 
Be pe’ (ab) ¥ (ib) +-pr2K'(b) 
The condition that the velocity vector is parallel to the magnetic induc 
tion vector Is 


| ¥ (hs) 4-(1 4or)pr?y"(b) A (yp) 
pe’ pe’ (yb) ¥ (us) 4 pr= hk’ (yp) . 


; _ 
so that 2K (pb)! 1- YT ie 0, 
f | 477 | 


or pr2K' (yb). 0. 
Accordingly the velocity and magnetic induction vectors will be parallel 
when (a) Kb) =0, 

(b) A 0. 
The first case is that in which the electric intensity field E is identically 
zero; if also the meridional flow vanishes, leaving only a rotational flow 


and a transverse magnetic field, then from (1.23) 


porp’(0) : r Bo. 
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This equation relates the angular momentum of a gaseous ring to the work 
done in conveying a unit pole around the ring and is an expression of 
Ferraro’s law of isorotation. Bernoulli's equation in this case is affected by 
the presence of the magnetic field through the term 4v?. The second case 
in which the velocity vector is parallel to the magnetic induction vector is 
when ./ = 0, and this case corresponds to the appearance of a singularity 
in the equation (1.36). When .o/ = 0 equations (1.21) and (1.22) will be 
consistent only if “ KW) 

- er Why’ 


2 A'(p) 
¥(p) 


Hence y is a function of r? only and the meridional flow must be entirely 


that is, pe’ (pb) = —r 


in the z-direction, hence the complete flow will be a spiralling motion on 
the surface of concentric cylinders r = constant. 

When ,’ is a constant (u’ = po) and © is a constant, equation (1.36) 
becomes 


1 {é sis ae (= =) = Fp) =m, “y ‘(), (1.37) 


pr\cz\pr ez) * ér\pr er]| 


and in this case we can interpret the motion as a rotational flow of a com- 
pressible fluid whose effective density P is modified from the original 
density p by the presence of the magnetic field, the relation between P 
and p being 

5 l af ‘ 
=—. (1.38) 
Pp 


If we use the earlier definition of ./ we have 


ae: La ee Se (1.39) 
1—(147r)pp'? 1 —p/py 
where uj; — 47/py, and thus the effective density will be increased when 
0 <p <p, At the lower end of this range, which will occur in a gas 
which is near the state of cavitation, the effective density P and true 
density p are closely the same. When we have a fluid whose density p 
is near p, it would appear that the effective density can be very much 
greater than the true density, the fluid density therefore becomes aug- 
mented due to the presence of the magnetic field. 
When yp’ is not a constant the above interpretation is no longer possible 
owing to the presence of the terms p’«” on the right-hand side of (1.36) 
which cannot be readily combined with the terms on the left-hand side. 


5092 .50 N 
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3. Case of incompressibility 
It is important to point out that equations (1.32) and (1.36) in the in- 
compressible case p = constant = p, can be written in the form 
y § k'( al 
P+ 4g +9) By = Fy), 
Pe 4a 
] 2 ele “(= 2) 
r cr; C2 | 


9 
p-r\er r dz 


F'(b) — Kk’ (yp) rBp _v ”"(h) +n") wee P (wet W?), 
4n r T T™p,.T 


1 
where Af = 1 p.m *(p). 
4n 


Accordingly the singularity ./ = 0 is equally present in incompressible 
flow. A certain amount of simplification can be effected in the % equation 
in this case. We have + # 

oe on 

cr 27 «Or 
and similarly for ¢.o//éz. Hence 


L (é pen ts Poh nt Ot OF Oa 
p2rier\ r cz cz\ r éz]} 2p? r?| er or" éz G2) 
” K" (pb) v , r Bg ” 
r ines Wet 

(yb) as rBg “( (ys) ie B (p) 


The left-hand side of this equation can be written in the form 
Ail ¢ (=> 1 (<2 | 
p> | Yr Cz 3 


and thus the complete equation could be interpreted as a compressible 


r oz ror\ r Or 
flow in which the density is proportional to </~*. More conveniently, 
however, we can introduce a new stream function ¥ related to the original 


one % by the equation 
4 
Y= | Wi dy 


I 9 , 
-ptu'™(yy\ tay, 
4a 
and in terms of this new function we have 
Oyre(lev\  te(ver ’ k" Uly- , "Bo » 
ot | c (- c \4 c | AS — ff (ys) ——— r Bo —<9 -}. —! as ()). 
p2r\r éz\r ez rer\r er |] 4a r 4r 
Accordingly the actual incompressible flow can be considered to be related 


to a second irrotational incompressible flow in which the stream lines are 
differently numbered. 
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4. Simple compressible flow solutions 
Case 1. Suppose that 
V = ui(r)k, B= B,6, j = jek, p=plr), p' (pb) = 9, 
Fis) = Fy, Ol) =O, K'(b) = Ky; Ky, Oo. Ko are constants. 
In this case we have rBy = pr*Ko, 
] 
w= — dep 
pr dr 
since the above conditions imply that V(y%) = 0. Hence we have from 
(1.31) ; pane 
Yo oy 14 Jw? + prKo = F,, 
-] ] 4n 
and from (1.36) it follows that w is constant and equal, say, tow,. If p has 
the value p, on r = 0 the above Bernoulli equation gives 
Yo py 4 pr?Ky a Yo px 1 
y—l1 4a y—1 - 
It is quite easy to show in this case that p, is the maximum density and 
that p tends monotonically to zero as r increases to infinity. As r-> 


p~ 47 Pe ‘ 

(y—1)kKG or? 
and it then follows that Bg is zero at r = 0 and tends to zero as r-> @, 
such that 


we have 


>. 

(y—1)Ky +r 
Hence By, attains a maximum for an intermediate value of r. The current 
j. is given by <a 


"Iz An sr Ko): 


j. is a maximum on r = 0 and tends monotonically to zero as r -> 00. 
Case 2. Suppose that 
V = vir)b4 w(r)k, j - j.(ryk, p=pr) B= B,(r)6, 
pw =0, Fb) =F, Ob) =O, Kb) = Ky V(b) = %. 

This is the same as in Case 1 with the addition of a rotational velocity, 

and we have from (1.23) and (1.24) 
wr=7,, 

Bor = pr*Ky. 


Bernoulli’s equation (1.31) becomes 


4 -) 1 . 
¥ 7 p’1+-4(v?4 w*) +7 Kor Bo _ 
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and the vorticity equation (1.33) or (1.36) becomes, as in Case 1, 


1 dp 


j constant. 
pr ar 


Since w and F, are constants, Bernoulli’s equation can be written in the 
_ ei, 74, Fh 
pyY-* + —— + pr constant. 
l 2r* 4 


Y 


It is clear that in this case, if there is no rigid internal boundary a vacuum 
will develop around the axis r = 0, and that if the radius of this vacated 
region is r9, then the value of the constant on the right-hand side is 


¥ 4% 2r5. Hence, as r >» %, the appropriate asymptotic expansion for p 


will be mC oe 
ony oT 


kz 


pr 


With p vanishing at r — r, and r © it is clear that it will attain a maxi 
mum value at some r between 7, and «. This constitutes the magneto- 


hydrodynamic circular vortex. 


Case 3. Suppose that 

V = v(r)6 4 wir)k, j= jo(r)6+ jk, p= pir), B= B,(r)6+ Birk, 
ise tr . , - ; . 
x constant, F KF, 9 ©,, AK’ (pb) 7 0. 
Po 

This case is a different type of extension of Case |: whilst there is a 
rotational component of velocity, as in Case 2, its origin is here connected 
with the non-vanishing of «’ and K’, that is, it is connected with the 
magnetic induction field rather than with the shear velocity case in which 
¥ is non-vanishing. 

We thus have from (1.23) and (1.24) 


| : 
a u'pr*K,, 
| = FT f 0 


(1 ,) Ba pr? Ky. 


It is interesting to note that in this case 


B. By 


v 


hence 





ON STEADY AXIALLY SYMMETRIC SOLUTIONS 181 


We shall assume now that B, and wu are identically zero and, with the 
above conditions upon the parameters, we obtain from (1.36) 


“(1 —_ e\n = constant. 
pr Po/ @r 


x ld 
Since dy 


= w this may be written in the alternative form 
pr ar 


w(1 a - constant = (1 -*, 
Po Po 


where w,, p, refer to values on r = 0. 
Bernoulli's equation (1.31) now becomes 


-) 
Oa ptt det Jet «Kur By = Fy 


Substituting in this equation for w, v, and Bg it becomes 


~(P%/Po)}* pr ks _ __ pre Ks 


(p/po)}® ' Spo —(p po} | 4x1 —(p/po)} 


which simplifies to 


Yo py-t4 weil — (ps Po)}* L pr? K G2 -(p Po)} 
a 2{1—(p/po)}? 8 1—(p/py)}? 
This relation defines p in terms of r and is the appropriate generalization 
of the corresponding result in Case 1 which can be deduced by making 
Po > &. 

Using conditions on r = 0 the value of F, becomes 


“) 
F, Y% px 14 hwy, 
yl 


and hence 


rk? (yop vs. ths aes 9) Px\” 
; a . oe ( Y Y 14 ue ( a ae 31 3) | 
Sir pi2—(p/po)}lly—1 a r ™" Po ws Po 


Various cases arise as follows: 


(i) O< px < po. 

The interesting feature of this case is that the solution for the density 
in terms of r becomes non-unique for r sufficiently large, in fact for r > rp, 
in the diagram. The interpretation of non-uniqueness is not to be obtained 
from the present theory but it may be connected, as in other examples in 
compressible flow, with the appearance of shock waves in the real 
problem. 
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i ; a 
(11) po < Pa < 2po- 
The unique solution in this case occurs only when r, <r < rz that is 
within an annulus. 





| 
{ . | 
_ _ - —————— H —_ —_ —s - 
lic. 1. Graph of p (abscissa) against Fic. 2. Graph of p (abscissa) against 
r?K2/87. Case i: 0 < pe < po- r?K2/87. Case il: py < pa < 2po- 





L 4 —— 


3. Graph of p (abscissa) against r?K7/87. Case ili: py, > 2po- 
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(ili) py > 2pp. 

In this case there is no unique solution to the problem. 

Accordingly when a singular point is present in the solution to this 
problem, a steady solution is only possible under certain circumstances, 
namely, when the axial density p, is less than some fundamental density 
2p) whose value is not determined by the theory. In the case, 0 < py < po 
there is no infinity in the velocity field since p = p, lies outside the domain 
AB and the solution furthermore is limited to a cylinder 0 <r < rg. 
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SUMMARY 

A study is made of the unsteady flow generated in a viscous heat-conducting gas 
by an infinite flat plate set suddenly in motion in its own plane. The corresponding 
Navier-Stokes equations, linearized for small Mach number, are treated by opera- 
tional methods and yield (1) a large-time representation of the flow, (2) a short-time 
representation, (3) a method for numerical evaluation of the flow at arbitrary time. 
The solution exhibits a transition of the flow from its short-time to its large-time 
form, At very short times the solution coincides with the initial-motion approxima- 
tion given previously by Howarth, and is characterized by sudden initial rises in 
pressure and temperature, At large times the flow splits asymptotically into two 
distinct motions, one of which is related closely to the incompressible solution, while 
the other consists of an acoustic wave modified by viscosity and heat conduction. 


1. Introduction 

{AYLEIGH’s problem in fluid dynamics relates to the evaluation of the 
unsteady flow which arises when an infinite flat plate, submerged in a 
viscous, originally quiescent fluid, is impulsively started moving in its 
own plane with constant velocity. This flow, though simple for an in- 
compressible fluid, becomes complicated when compressibility is taken into 
account. The heat generated by viscous dissipation produces in a com- 
pressible fluid not only temperature variation but also variations of 
pressure and density, and thus induces a velocity component normal to 


the plate. The corresponding Navier-Stokes equations are highly non- 


linear and cannot be solved exactly. 

To attack the problem for a viscous gas, Illingworth (1) and later van 
Dyke (2) used the boundary layer approximations. A more unified ap- 
proach was initiated (somewhat earlier) by Lagerstrom, Cole, and Trilling 
(3), who linearized the equations of motion under the assumption that 
the Mach number M is small. However, the solution given by these 
authors neglected the thermal conductivity of the gas, and was expressed 
in terms of repeated integrals in complex form which were not evaluated. 
Further progress was made by Howarth (4) who took thermal conduc- 
tivity into account and, applying operational calculus to the linearized 


[Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 2, 1960} 
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equations, obtained the Laplace transform of the pressure. This trans- 
form turned out to be very involved, but Howarth succeeded in deducing 
from it the variation with time of the pressure at the plate. Howarth then 
returned to the non-linear equations and derived an approximate solution 
describing the beginning of the flow. Apart from this initial solution the 
flow in the interior of the fluid remained unknown. 

The purpose of the present work is to complete the solution of the 
linearized Rayleigh’s problem for a heat-conducting viscous gas. As in 
Howarth’s work, the Prandtl number is assumed to be } and the plate is 
assumed to be heat insulated. The linearized equations of motion (which 
represent the first step of expanding the full solution in powers of M*) are 
treated by the methods of operational calculus. Laplace transforms are 
obtained of the pressure, temperature, and density variations and of the 
normal component of velocity. These transforms are first simplified, and 
then interpreted to give the following results: 

(I) an asymptotic representation of the flow, valid at large time (and 
arbitrary distance from the plate); 

(Il) a short-time approximation, which coincides with Howarth’s initial- 
motion solution; 

(111) a Fourier-integral representation of the solution, which facilitates 
the numerical evaluation of the flow at any time and distance. 


A discussion and a graphical illustration of the flow are given in the 


concluding section. 


2. Equations of motion 
The quantities to be evaluated are the pressure P, density p, tempera- 
ture 7’, the velocity u parallel to the motion of the plate, and the velocity 
v normal to the plate. As the plate is plane and infinite, these quantities 
depend only on the normal distance y from the plate and on the time t. 
The corresponding equations of motion for a viscous, heat conducting, 

perfect gas are 
eu é[ eu 


ou . fe 
Plat ey) oy\" by p 


Dee 2, .p ‘ Des 
(B+0%) — S48 F(.) (1b) 
ot cy cy 3 ey\ cy 


(la) 


= Td yy = OC = le,d 
PS Gy (PY) , P= RT, (le,d) 


PF . 2%) oP _ oP 2/28). [sy Sf\"], (ney 
. Ca ey) et ty by ey] TP Ney) * 3hey) 
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where ,: and « denote the coefficients of viscosity and of heat conduction, 
R is the gas constant, and (, the specific heat. The quantities » and « 
depend only on the temperature, and x = (C,/Pr)z where Pr denotes the 
Prandtl number. At time ¢ = 0 it is supposed that the thermally insulated 
plate starts moving in its own plane with the uniform velocity U. Indicat- 
ing by the subscript 0 the undisturbed state of the gas, we have the 
following boundary conditions: 
when ¢ < 0 (for all y), and as y > » (for t > 0): 
ne 
Fi T) Po 


at y © (for t > 0): 


, T 
U, : 0. y= 0). (2b) 
oy 
The equations of motion and the boundary conditions permit us to expand 
the solutions in powers of M?, where M = U’/a, and a, denotes the inviscid 
speed of sound. When the Mach number is small the solution has the 
form 
= ay| Mu,+O(M?*)], v = a,| M*v,+-O(.M‘)], 
P = Pi[l+ M*P,+0(M‘)), T = T,{1+ M?T,+ O(M‘)], 
P Poll + M*p, T O(M*)). 
For u, the equations of motion give 
cu ctu 
= v%—~ = 0, 


ao 
cy* 


where v) = 49 Pp denotes the kinematic viscosity. Therefore 


? 
‘, = erfe| J ) (5) 
» : t) 
2, (Ho 
which is identical with the velocity solution in the case of an incompressible 
fluid. Furthermore, we assume that Pr = ? and we introduce the dimen- 
sionless distance and time variables 
3 a? 
: —é, 


(6) 
4y M% 


where y denotes the ratio of the specific heats. From the equations of 


motion it follows then that r,, ?,, T,, and p, must satisfy the linear system 
of equations 


ev, @P, ev - 
y— i += 0, (7 a) 


’ Of on cn 


0, P—T,-p, =, (7 b, ¢) 
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eT, CP, OT,  Vyv— D1, ey sxnen (7d) 
tor c 2 


€ Cr CH 7 T 


¥ (y—1) 
3. Operational solutions 

The form of the boundary conditions makes the linearized problem 
suitable for treatment by the operational calculus. Denoting by f(7, p) 
the p-multiplied Laplace transform of f(y, 7), we apply the Laplace trans- 
formation to equations (7) and obtain 


vp, + Ms _. 0, (sa) 


PP. +5, = 0. P,—T,—p, = . (Sb, ¢) 


a2 
ypT,—(y— 1)pP. = = ue = bpk,(anvp) (8d) 
. 


where A, denotes the modified Bessel function, and where 
2 Sy! 
b- y(y—1), = (7) 
7 3 


The boundary conditions become accordingly 
rs =_ P, -_- T, = Py —- as 7) >, (10a) 
... 
i,.= - 
oy 
To find P, we eliminate 7,, f,, and *, from equations (8), getting 
eP, Pp 


A : sf 
ey? a itp! = ere (11) 


=0 atyn= 0. (10 b) 


The relative simplicity of this differential equation is due to the choice 
Pr = }; for any other value of Pr the resulting equation would contain 
the fourth derivative of P, with respect to 7. The boundary conditions 


for equation (11) are 
P,=0 asn>, 2-0 aty=9, (12) 
cn 
the second of these requirements being a consequence of the conditions 
(10b) and the equations (8a,b,c). The solution of equation (11) by the 
method of variable parameters yields the Laplace transform of the pressure 
variation: 
7 a 
P, = bq E a7 i) cosh(qA)Ko(avp A) dA+-cosh(qn) ( e-WK,(avpa) aa| (13) 


0 7 


P (14) 


where ee aa 
. v(l+p) 
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This operational form was first obtained by Howarth (4) and was inverted 
by him for 7 = 0 to give the pressure at the plate. Inversion of this 
transform for 7 > 0 was, however, not undertaken. 

The transform of v, is found by solving the equation 


a , 


2 2 


(yp+q)— 
ype 
7* YET oy? 


Lypq?t, = tet K,(anvp), (15) 


which is the result of eliminating p,, 7,, P, from equations (8). The 
boundary conditions for *, follow from the conditions (10) and the equa- 
tions (Sa), (Sb), giving 


0 asn>D, .= — QO at y= 0. (16) 


The application of the method of variable parameters and a subsequent 
integration by parts lead to the operational solution for the normal velocity 
component: 


sje-™ ( cosh(qA) Kg(av pA) dA—sinh(qy) [ AK (avp aA) dA— 


0 U) 


q? 
yPp— 


e~“(ypm | cosh(.(yp)A)Ko(av pA) dA+ sinh(s(yp)n) [ e-“YPAK (avpa) dal- 


0 


7 = 
(17) 
The corresponding transforms of the density and temperature can now 
be determined directly from the solutions (13) and (17), since 
1 ©, P 1 oF, 
Ps ,, A Pp, +-— ~ S 
P € 1) Pp cn 


but there is no need to state the resulting expressions explicitly. 


4. A modified form of the operational solutions 

The operational solutions obtained above cannot be readily interpreted 
in their present forms, but it is possible to derive useful results from them 
after certain simplifications have been introduced. One such simplification 
is achieved by substituting the formula 

A,y(as pa) ( e-a%pAcoshy dy (18) 
0 

((5) 181) into the transforms (13) and (17), and by subsequently changing 
the order of integration. It is then found that the integration with respect 
to A involves exponential functions only and can be readily carried out. 
In some of the resulting terms the integration over v can be also per- 
formed, and we get 


D> p p ae aes ~ ( 
P, 21 7 492; Va1 + Vat Veg, (19) 
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where 
P = zal Pp : | fem | 


1 = 3 att (ae I)p 


D ’ 
e-anvpcoshy 
P,. = —bp [ dy 
J a®(1+p)cosh*y— p 
0 


m™ Be F meme 





2 y+ (y—Vpla?+ (a?—1)p 


mh ‘Ps g-gvom 


(ay) y+ (y— Dp 


baxp 


x 
[ e-anvpcoshy cosh v dv 
| ®(1 + p)cosh*y— p |(a? cosh*y — ”) 





0 
We now observe that the transforms P,, and @,, are of the same type as 
those investigated by Lagerstrom, Cole, and Trilling (3), and by the 
present author (6), in connexion with the propagation of longitudinal 
disturbances in viscous, compressible media. On the other hand, #,, and 
the integrands of P,, and *,3 are operational forms of the kind encountered 
in problems of heat transfer (or, in fact, in Rayleigh’s problem for incom- 
pressible fluids). Inversion of the operational solutions thus appears within 
our reach. This remark applies not only to P, and @, but also to 7, and pg, 
for we have 
Po = Par > Poet Pas: T, - Ty, Ty, 7 Ts, 
where, as follows from equations (19), (21), and (8b), 


h 


mb l Pp 2, wisirsp> | 


2 70TH Te 
‘ a... 2 : te — e-n(yp) 
2y(0?—y) y+(y—Dp 


Pre 


x 
- 


. bx? | e-anvpeeshy cosh*y dy 
P23 ax ee a ke ee 
| a?(1-+ p)cosh*v— p |(a®cosh*y —y) 





and where the corresponding terms of 7, are given by 7, = P,—A,. 


5. Asymptotic solution for large time 


Starting from the operational solutions in their modified forms, we pro- 
ceed to derive asymptotic representations of the actual solutions, valid 
for large values of the time variable 7. We consider first the transform 
P,,. Its interpretation is found, according to Mellin’s theorem, by 
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evaluating the complex integral 
6 eTAg—nA vl +ay x2 + (a? 1)A] 4 dA (24) 
= J \A 


cia 


P,,(, T) 


Following an analysis given in (3), we deform the path of integration to 
the contour composed of the unit circle with centre at A = —1 and of 
the two straight lines stretched between A = —oo and A 2 just below 
and just above the negative real axis. It appears then that for large + 
the only important contribution to P,, comes from the are —6, < 6 < 4, 
of the unit circle A = —1-+e'®, where 6, is an arbitrary positive con- 
stant smaller than 7. In fact, the contributions e, of all the other parts 
of the contour are found to vanish exponentially when 7 becomes large. 
Thus, using s = sin(4@) as the variable of integration along the arc, we have 
8; 
Py 2-8 * 278? +2i(r—napith(s)a/ 9) wes L€, 


: v(is) | 


8) 


where 0 < s, < 1, and where 


h(s) = 2s{,(1—s?)—1} ) (26) 

g(s) = [a2+ 2i(a®— 1)s{(1—s2) is} }-' (1 —s?) + is}t(1 —a2)-4 | 
The functions g(s) and h(s) are now replaced by their values at s = 0, and 
the resulting errors are estimated. When s) < 8, both g(s) and h(s) are 
analytic, and in addition A(0) = h’(0) = h’(0)= 0. It follows that 
g(s)—g(0) < A's and that h(s) < Bs*\, where A and B are positive 
constants. Since h(s) is real in the interval —s, < s < s, we have further- 
more \¢'7®)—1) <7 h(s) < Br's*\. Therefore, denoting by «, the error 
when g(0) x! is substituted for g(s) in (25), and denoting by ¢, the 
error made when in the resulting integral e'7" is replaced by 1, we obtain 


$4 


ds 


; + €) + €g 7 €3- 
y (ts) 


where 
2-ibA | e- 278" g . 2 T(3)bAr : 


. 
"1 


ae 
€; <2 '~ Br ars"! g' =< 2 "1° Br i 


x 


now extend the limits of the integral in (27) from ( 81, 
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(—o0, +00); the error caused by this extension vanishes exponentially at 
large r. It follows that 


b [ e-sntsste "8 ds or-4) 
y (ts) © 


zx 
a 


2 , ds 
e *7*l cos 2(7 — n)s + sin 2(7— )s| ee O(r-*). 
\8 
0 


To evaluate the last integral, we observe that 


[cos 2(7— y)s + sin 2(7—y)s] a v(7\n—7|)[F_\(2|n—7!8) FJ, (2|n—7/8)] 
v8 
(30) 
where + corresponds to 7 2 7. We recall then the formula ((5) 394) 


‘ae 


Vr .. 
e- C8 J (as) ds = —e-a'8f, | (31) 
! * 2c ty se? 


in which J, and J, denote the Bessel functions, p > —1, ¢ > 0,4 > 0. 
Combining equations (29), (30), and (31) we find that 


Pa(n.7) = Ter F() +O(r-) (32) 


where, as in the following, F denotes the function 
F(x) = e (a?) T(x?) FU,(2?)] (FF when x © 9). (33) 

The function F(x) is positive and bounded, has a single maximum point, 
vanishes like 2! when x —o and exponentially when x +. 
A graph of F(x) is shown in Fig. 1. 

We turn now to the interpretation of P,, for large +. Applying Mellin’s 
inversion theorem and subsequently changing the order of integration 
we obtain 


} x “+4u 


: dy . l : | eta aqvAcoshof j 1 Af(v)| 1 dy, (34) 
v | cosh*®y 277 


0 c—ta« 
where fe) = 1 -- ey (35) 
x? cosh*y 
The complex integral with respect to A can be taken along the contour 
consisting of the lower and the upper sides of the negative real axis (in- 
dented at A -1/f(v)). We find then that the only contribution to P,, 
which does not vanish exponentially for large + arises from the part of 
the contour near the origin A = 0. In the neighbourhood of the origin 


| 1+-Af(v)|-! can be approximated by 1, since | f(v)| < 1 forany v. Evaluating 
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the resulting complex integral, and estimating the error due to the 
approximation, we get 
P,o(n,7) = —- b rol L | ¢-atnttrcoshty ae O(r* 


2Zavr wry cosh v 


0 


bxrn 
. 7} " erfe( 2” 


4x NT 


\} + O(r-). (36) 
NT 
‘The required representation of the pressure solution P, is now obtained 
at once by adding P,, and P,,. Equation (36) shows, however, that Py», 
being of the order r~!, should be neglected at large times in comparison 
with the error term of P,,.. Therefore representation of the pressure for 
large 7 is - ee 
P(»,7) = — F(? ); O(r~*). (37) 
4x \ 2N7 

The interpretation of @,,, f,,, and 7;, for large + can be carried out in 
the same way as that of P,,, and yields representations proportional to 

‘F((m—7) 2,7) with error terms of the order of magnitude O(7~*). All 
the other transforms which appear in the modified operational solutions 
(21) and (23) can be interpreted either as for RP. or as for the integrand 
of P,,: the resulting terms cannot, however, be neglected in comparison 
with 7-'. Thus we find that at large + the normal component of velocity 
is given by 


7 . T x7 
Vo(9, T) a er - l _ 
y ‘ 2ay \2vT 


rh river exfe( = — 1M LO(r-2) (38) 
| INF 


} he 

2y,(a*—y) 
and that the corresponding representations of the density and the tem- 
perature are 


by7 
Pol, T) 


\ Ly (a* 


/ 


ry byr NY 7 TT 
T,(y,7) oN he ] e. 
Viyla*—y)} iim 


where G(y) = ¢ foal a) ds. (41) 


The function G(y) is introduced into these results by the integration with 
respect to v required for the evaluation of p,, and 7,,. Fig. 2 shows a 
graph of G(y) for y 1-40. 
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F (x) G(y) 











Fic. | 


6. Short-time approximation 

The modified form of the operational solutions enables us to determine 
also the behaviour of the flow at the beginning of the motion. For this 
purpose p may be regarded as large, and we notice that for large p the 
operational solutions of section 4 reduce to well-known Laplace transforms 
(thus e-7”*"+”) becomes e-7”). Interpreting these transforms by means 
of tables and integrating where necessary with respect to v, we obtain an 
approximate solution of the problem, valid for small values of +. We find 
then that the solution for small 7 coincides with the initial-motion solution 
given by Howarth (4). 

It should be pointed out that in Howarth’s work the initial-motion 
solution was derived from the full (non-linear) equations of motion simpli- 
fied by assuming constant viscosity, while in the present work the linearized 
equations of motion are employed. Still, the coincidence of the results 
is not fortuitous. In fact, when 7 is small the normal velocity and the 
variation of density are negligible; thus, if constant viscosity is assumed, 
the full equations of motion reduce for small 7 to the linearized equations. 


7. Representation of solution by Fourier integrals 

To complete the solution, a method must be found for evaluating the 
flow at unrestricted values of time and distance. As both the original 
operational solutions and their modified forms given in section 4 appear 
incapable of yielding such a method, we transform the operational solu- 
tions in yet another way. We substitute Basset’s formula (5, 172) 


V(a2p+-o?) 


K,(avpaA) = | eee ad (42) 
0 


into the original operators (13), (17); then changing the order of integra- 
tion we first carry out the elementary integration with respect to A. This 


procedure results in transforming the solutions into Fourier integrals. In 
5092 .50 O 
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fact the operational solution P, becomes 


P,(y, p) b ( d,(o, p)cos(ya) do, 


p* 


Ww here d io, P) = , = = —— mae 
. (p* + o*p+oa*), (a°p4 o*) 


it follows then that the actual pressure solution is given by 
P,(y,7) b | d,(o,7)cos(ye) do, (45) 
0 
é,(o,7) denoting the inverse Laplace transform of ¢,(¢,p). Similar 
Fourier-integral representations are found for the remaining solutions. 
Thus 
r z 
b | d,(o, 7)sin( ya) do, p(7,7) = 6 | $,(0,7)cos(na) do, (46) 
6 d 
», and d, denote the inverse Laplace transforms of 
op" Co 
: ; di < ———s d (o, p) 4,0, p). 
(yp o*)( p* rop- a~)x ( “p+ a”) " p 
(47) 
The corresponding result for the temperature variation is obtained, of 
course, by using 7, P,—p,. A glance at equations (44), (47) shows now 
that the operational forms d(a, p) are sufficiently simple to permit in- 
terpretation. This provides a means for computing the actual solutions 
numerically at any given values of » and r. To compute the solutions we 
evaluate first the functions 4(¢,7) by interpreting the Laplace transforms 
é(o,p), and then we calculate numerically the Fourier integrals (45) and 
(46). 
The interpretation of d(o, p), though elementary, is somewhat tedious. 
Applying to 4, the rules of operational calculus and writing 


21 ,() lala ; , (a" $)}, ry 9(e) fq 2g? 21 9(0)} (48) 


we get 
2 eT erf(r, v7) |. (49) 
ry re 


> 


| ee 
po, 7) i “7 erf(r, v7) 
XA< — <9) | 


When o > 2 the quantities 7,, r, are either real or pure imaginary, and 
then ¢, can be computed from tabulated functions. However, when o < 2 
equation (49) ceases to be useful, for then r, and r, become complex; in 
this case it is best to evaluate 4, with the aid of suitable expansions. 
A power series for ¢,, useful for small values of o?r, can be deduced 
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directly from the transform (44). For large values of or an asymptotic 
expansion of ¢, can be obtained from (49) by employing the well-known 
asymptotic expansion of the error function. Evaluation of the other ¢ 
functions proceeds in the same way, the additional factor in (47) intro- 
ducing no essential difficulties. All the functions ¢ vanish quickly with 


P, $ 


-- 











Fia. 4 
increasing o (and 7), so that computation of the Fourier integrals is carried 
out readily. 
8. Discussion and numerical results 
The solutions P,(,7), V2(7,7), p2(n, 7), computed for + = }, 1, 4, and 10 
by using the Fourier-integral representations, are shown in full lines in 
Figs. 3-5 and the solution 7;,(y,7) computed similarly for r = 1, 10, 25 is 
shown in Fig. 6. The dashed lines show the asymptotic solution obtained 
in section 5 and indicate that the large-time representations provide 
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a good approximation when 7 > 4. The initial-motion solution is found, 
on the other hand, to be valid for > << } only. For all the values of 7 the 
pressure at the plate (» = ©) calculated from the Fourier integral agrees 
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20 30 ” 


Fic. 6 


very well with the values given by Howarth (4); this serves as a useful 
check for the computations. A complete agreement exists also between 
Howarth’s formula for the large-time variation of the pressure at the 
plate and the corresponding result obtained from the present asymptotic 


solution. In fact equations (37) and (33) give for the pressure at the plate 
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at large + 


mb mb ; 


P,(0, 1) © rt F(— 301) = [1 y(t7)+4,(47)]: 


«4 
2a 
replacing the Bessel functions by their asymptotic forms we get 


(0,7) 4, (51) 
~~. 
which coincides with Howarth’s result. 

The outstanding property of the solution is the transition of the initial 
motion into the large-time flow. In the initial stage the normal velocity 
and the variation of density are still very small; we have then P, = 7), 
so that the dissipation of the main motion raises the temperature and the 
pressure equally (cf. 4). The initial variations of the temperature and 
pressure resemble the temperature variation in the incompressible Ray- 
leigh’s problem: they decrease from a constant value at the plate to zero 
in a narrow but expanding region. During the transition period the 
character of the flow is changed by the appearance of a wave-like motion. 
The pressure variation is transformed into a progressing pulse, and a similar 
normal-velocity pulse builds up. A density disturbance is developed, con- 
sisting of a condensation wave moving forward and a rarefaction ‘tail’ 
ending at the plate. The temperature variation (Fig. 6) exhibits also the 
formation of a pulse ahead of the expanding boundary region. 

An insight into the large-time flow pattern is obtained by examining 
the asymptotic solutions of section 5. These representations have the 
forms 


P, ~— PH; v~ vf) 4 vi); p» ~ py!’ + py”: T» sla T j T?, (52) 


where upper affix (1) denotes the terms containing F and affix (2) refers 
to the remaining terms. Each of these sets of terms has a definite physical 
meaning, so that the large-time flow can be regarded as a superposition 
of two distinct flows. The flow described by the terms with affix (2) is 
related closely to the incompressible solution of Rayleigh’s problem. Thus 
PY = 6: moreover, TY) coincides exactly with the temperature solution 
in the incompressible case, since it satisfies the corresponding energy 


equation orp 2p 2 
eT vy o?7 Vp [Cu — 
0 “( ) (53) 


at Prey ec, \ey 

and the boundary conditions. The density variation and the vertical 
velocity do not vanish however, for the equation of state requires that 
p,” = TY, from which the vertical velocity vy? is then determined by 


continuity. It may be noted also that the terms with affix (2), when 
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expressed by the physical time and distance, are found to be independent 
of y. 


The terms with affix (1) describe, on the other hand, the propagation 


of an aperiodic wave in a compressible fluid, modified by effects of viscosity 


and heat conduction. This wave is isentropic, for we have 


Py) = yp), TY = (y—1)p2? (54) 


furthermore, the terms with affix (1) asymptotically satisfy the equations 
of motion with the dissipation term omitted. Dissipation thus determines 
the wave indirectly by producing the generating disturbances at the plate. 
The wave (being derived from operators containing e~””*"*?)) is of the 
same general type as the motions studied in (3) and (6). The results of 
section 5 show that the wave produces positive variations of all the 
quantities P, 7’, p, and v, and that these variations form a single propagat- 
ing pulse. At large times the peak of this pulse moves with a speed 
approaching the speed of sound, the peak amplitude decreases like 7! 
and the thickness of the wave region increases like r?. 
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SUMMARY 

It is shown that the full equations of motion for the steady flow of an incom- 
pressible liquid have a siraple approximate solution, valid in a boundary layer near 
the plate y* < 4ar, z= 0, when the velocity at infinite distance from the plate is 
of constant magnitude U and parallel to the a-axis, This solution, which is fully 
three-dimensional, is compared with Moore’s approximate boundary-layer solution 
(1); and on the basis of the comparison suggestions are made about the validity of 
Moore’s solution in the case of flows past plates with more general shapes for the 
leading edge. 


1. Introduction 


Tue well-known Blasius boundary-layer solution is not the only form of 


boundary-layer solution to the two-dimensional problem of flow past the 
plate x > 0, 2 == 0, when the velocity at infinity is of constant magnitude 
and parallel to the x-axis; for Lewis and Carrier (2) have shown that 
another boundary-layer solution is obtained by transforming the equa- 


9 


tions to parabolic cylinder coordinates €, ¢, defined by x = &—2@, 
2 2é¢. By using these coordinates, Stewartson (3) and others, whose 
results are quoted by Kaplun (4), have shown independently that the 
Oseen linearized equations for this two-dimensional flow have an exact 
solution. In a previous paper the present author (5) has also shown that 
the Oseen equations for the three-dimensional flow past the plate y* < 4az, 
z = 0, when the constant velocity at infinite distance from the plate is 
(1,0, 0), have an exact solution when they are transformed to a degenerate 
form of paraboloidal coordinates, defined by the equations 


x = a(€+7+ C—1), 
= 4a®(E€—1)(1—n)(1—Q), 
= —4a*én@, 
where & ‘3 = Nn = 0,0 = i. 
The boundary-layer solution of Lewis and Carrier and both the two- 
and three-dimensional solutions of the Oseen equations have one feature 


(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 2, 1960] 
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in common: the cartesian components of the velocity are of the form 


(x(Z), 0, 0) +-grad 4(Z), (2) 


where the scalar functions y(¢) and 4(¢) are functions of a single curvilinear 
coordinate ( in each case and have different forms in the three solutions. 
As the three-dimensional problem is merely a generalization of the two- 
dimensional problem, and the coordinates of equations (1) may be regarded 
as a generalization of parabolic-cylinder coordinates, it is natural to in- 
quire whether the three-dimensional problem of the plate with a parabolic 
leading edge has a boundary-layer as well as an Oseen solution of the 
form (2), when C is defined as in (1). 

The coordinate system defined by (1) is orthogonal and, if the length 
element is given by 

ds? h? d&* 4 hedy? + hdl, 

then 


p2 a*(€—n)(E—-C) 2 a*(n—C)(E 7) h2 — a*(€—L)(n—C) 
i 


(3) 
All the parametric surfaces are paraboloids with their axes along the 
x-axis. Those for which € — constant have their vertices pointing in the 
direction of the positive x-axis and are cut by the planes x — constant 
in ellipses, while the surfaces » — constant, which are more complicated 
in shape, are cut by the planes x — constant in hyperbolas. The surfaces 
¢ — constant are of some importance. If ¢ = —y (y constant > 0) the 
surface is the paraboloid 


4a(x + ay) 


which has its vertex pointing in the direction of the negative x-axis and 
is cut by the planes « — constant (> —ay) in ellipses with major and 
minor axes parallel to the y- and z-axes respectively. The paraboloid sur- 
rounds the plane area z = 0, y?—4ax < 0, which is the parametric surface 
{= 0. In fact the paraboloid encloses a layer round the plate past which 
the flow is taking place, the layer progressively increasing in thickness as 
y increases. The area in the plane z — 0 outside the plate (i.e. the area 
2 O 4? —4dar () is the surface » 0. 

If a point, whose cartesian coordinates are (x, y,z) and whose curvilinear 
coordinates are (€;», (). is moved so that its distance from the plate tends 
to infinity then ¢ - —#, and conversely: but if €-» —o with € and 7 
fixed then x >» —x. In order to increase the distance of the point from 
the plate to infinity in any other direction than that of the negative x-axis 
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it is necessary to impose a relationship between ( and é as ( -» —o. For 
example, if the point moves on the z-axis and its distance from the plate 
tends to infinity then 7 = 1, €+€ = 0, and [+ —m. 


2. The equations and their solution 

It may easily be verified that a vector with cartesian components 
(A,0,0) has components (a4 /h,,aA/h,,aA/hy) referred to the coordinate 
system (1). Thus the assumption of the form (2) for the velocity com- 
ponents in cartesian coordinates is equivalent to the assumption of the 


form Vaf() Vaf() Val f(¢) +-9(0)) 
2h,” 2h, | 2hs 


(4) 


for the components of velocity referred to the paraboloidal coordinates, 


with ‘ ‘ 
f() = 2x(¢) and = g(f) = : a 
l Va dl 

If these components of velocity are to satisfy the conditions at infinite 
distance from the plate (where the cartesian components of velocity are 
(U, 0, 0)) then the components given by (4) must tend to U'a/h,, Ua/h,, Vash, 
as (-» —%, in whatever direction distance from the plate is increased, 
This condition will be satisfied only if f({) -> 2, g(0) > 0 as €-+ —m. (The 
sufficiency of the condition is clear. Its necessity may be seen by letting 
distance from the plate tend to infinity along the z-axis.) 

The velocity components must all tend to zero as the plate is approached 
from any direction, and as the plate coincides with the surface ( — 0 this 
condition requires f(0) = g(0) = 0. 

The complete boundary conditions for f() and g(Z) are thus 

f(0) = gO) = 0 and f(l) +2, 9(f) >0 as C+ —m. (5) 

If the velocity assumed in (4) is to satisfy the equation of continuity 
divu — 0, then 
d 


j 
dt’ 


When the equation of motion 


CUO IO) + 00} FO) FILO = 0.6) 


F l 
u « curlu— 4 gradu? grad p+vcurlcurlu 
p 


is transformed to the €, 4, ¢ coordinates, and the assumed form for the 
velocity substituted, the equations which arise are 


U2q2 ri f+9) 4 U2a*g(2f + q) © ] ' 1 Op | vua “( h, P 
th, h2 . sh, o€\h3 ph, o€  2hyhy ol\hgh, ); 


(7) 
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1 ep via ¢ 


*a*g(2f+ 
Ais) 2h. sh Akt’ 


xf p a 7) 
, UFa*g(2f—g) ¢ (1 
+3 Sh, Aja) 
l cp rapa “2 | 
h, 


ph, ¢ C 2h, | icf 


ick) 


and g° denote df dt and dq df. The second term on the right- 


1 side of equation (7) may be transformed to 


wma 
é{1\) 
| ] . 
c h2 } 
3 


the right-hand side of equation (8) may be 


on 


d 
1—{)}* dc 


If now the press ire Is expressed is 
T2g25 vl ap pe 
j2 ayf—-g)— sia 
ths aN3 
)} may be written as 


th? cP 2y | 


3 
pl *a* c& la\ c(l 


ug eP =" | bg —eyy’y)._ (13) 
ta2@ cl Val{—G1—Zt)}' dl ° ies 

The left-hand side of the last equation is obtained by using equation 

6), while the right-hand side is obtained directly from equation (9) after 
some manipulation 

Were the velocity components (4) 

the functions f(f) and g(f), and the 

have to satisfy the four equations (6), (11), (12), 


an exact solution of the problem, 
pressure function P(£,,¢), would 
(13), the boundary condi- 


and P(é,»,¢) >a constant as [-» —o. In fact this system of 


However, for a high-speed flow it is to 


tions (5) 
equations has no exact solution. 
be expected that there will be a boundary layer near the plate, and that 
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all significant changes in velocity and pressure will occur in this layer. 
In that case f({) and g(Z) will undergo finite changes as ( varies from zero 
at the surface of the plate to —e, say, at the edge of the boundary layer. 
« being small. If the variable is changed from ( to ¢,, where ¢, le, 
it is clear that if f and g are O(1), then df d{, and dg d{, are also O(1), 
since the main changes in f and g occur as {, varies between 0 and —1. In 
this case df df and dg df are O(1 €), so that the term f‘(f+g) on the left 
of equations (11), (12), and (13) is O(1 €), while the terms g? (€—¢) and 


gq? (» —C) are O(1), provided 
g° (y—¢ ), J ne> I. (14) 


Thus, if the changes in the flow are confined to a narrow boundary 
layer near the plate and if the condition (14) is satisfied (as well as in 
certain other circumstances which will be described later), the second 
term on the left-hand side of each of the equations (11), (12), and (13) can 
be. ignored in comparison with the first. These equations can then be 


solved with P(£, yf) equal to a constant and with f(¢) and g(C) solutions 


of the simultaneous ordinary differential equations (6) and 


‘ _~ 
((f+9) = =~ Fl 
; 9g Va{[—i— Jide 


(1 —f)]*f’}. (15) 


It is convenient now to detine a Reynolds number R for the tlow, by the 
relation R lay, and to change the variable from ( to ¢ (= 0), where 
t R'}—U1—f)}'. The function F(t) is then detined as 

Fit) = RY —CIl—f)) fl) —9i0)). 
In terms of the new function and variable, equation (6) becomes 

f(Q) = dF dt, 

and equation (15) then simplifies to 

_,a2F d on OF 

F(t) —~ + itary? | 

dt®— dt| dt? | 

Since the velocity components (4) in the €, », ¢ directions can now be 
ladF ladF la 


2h, dt’ = 2h, dt” —hgt 


0. (17) 


written as 


Fit), (18) 


all the boundary conditions on the plate and at infinity will be satisfied if 

F(0) = (dF dt), =, and dF dt--2 as t+. (19) 

The condition for the flow to be of the boundary-layer type is that the 

teynolds number R should be sufficiently large. Now as R tends to 
infinity, equation (17) degenerates to 

fF @&F 


Fit . = 
(*) dt= ° dé 


(20) 
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which, with the boundary conditions (19), is the equation occurring in 
Blasius’s solution for flow past a flat plate with a straight leading edge. 
It is well known that the significant part of the solution of equation (20) 
is confined to the region 0 < ¢t < 4 (approximately); and so, for R large 
enough, the significant part of the solution of equation (17) will be con- 
fined to about the same values of t; that is, dF /dt will have approached 
very close to the value 2 when ¢ has reached a value of about 4. This 
means that significant variations in the velocity components (18) are con- 
fined to values of t less than about 4, or values of —C less than about 
16k". The flow given by the velocity components (18) is therefore 
certainly a boundary-layer flow for sufficiently large values of R, with 
a boundary-layer ‘thickness’ « (for f) which is O( R-*). It follows that the 
term 4R-'f? in equation (17) is O(e) in the boundary layer and can be 
ignored in comparison with unity, so that, according to boundary-layer 
theory, equation (17) can be replaced by equation (20). In fact the term 
containing R-! in equation (17) arises because {—({(1—¢)}! has been re- 
tained unaltered, whereas in the boundary layer —¢, which is O(e), is 
negligible in comparison with unity, although outside the boundary layer, 
for sufficiently large values of ¢|, the term 1—€ is clearly of a different 
order from unity. As ¢ > % in equation (17)—and one of the boundary 
conditions (19) is that dF dt + 2 as t + «o—the term 4R-'#? must be the 
dominant term under the radical sign. There is some interest in the reten- 
tion of this term, for a comparison of the solutions of equations (17) and 
(20) should throw some light on the way that one of the terms which 
would normally be neglected in boundary-layer theory affects the solu- 
tion. 

This comparison has been made for two moderate values of R. In order 
to avoid imposing a boundary condition at infinity as well as at the origin, 
equation (17) was transformed to 


, a*F. d | 
ea 
' ds? ds\ 


d*F,) = 


(1+ 4s? 
) ds? | 


(21) 


where F\(s) — R'F(t) ands = R-'t. The boundary conditions to be satis- 
fied by F\(s) are F, = dF, ds = 0 when s = 0 and dF,/ds -> 2R as 8s > o. 
By starting a numerical integration at s = 0, with F, and dF,/ds zero and 
d*F, ds? given an arbitrary value, the Reynolds number R correspond- 
ing to the chosen value of (d?F,/ds?),_, is found as }lim(dF,/ds). (By 


s--= 
choosing (d?F, ds*), , equal to 1-32824k!, a value for R close to k is 
obtained.) 
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Table showing values of F(t) and dF /dt found from 
equations (17) and (20) 





Blasius solution R = 99°4 8-49 
kiR=1 k = 100 9 





dk J dF 


i = F(t — 
“ dt sad di 





° ° ° 
0788 0239 7 orS15 
1458 0922 48 1°439 
1°347 1°g22 % 1-735 
1°976 3070 ‘97 I°g30 
1-998 4259 * 1-980 
2"000 5°454 , 5 1-995 
2°000 6651 “9995 1°gQ99 














2 
ti t-@ 














The integration has been performed for the two values k = 100 and 
k = 9. The corresponding values of R are 99-4 and 8-5. If the solutions of 
equation (17) for these values of R are compared with the solution of 
equation (20) it is found that (d?F /dé?),_, is increased by 0-9 per cent and 
9-1 per cent, respectively, and these figures also give the greatest per- 
centage change in both F(t) and dF /dt in the two cases. (The greatest 
percentage change in each case occurs as ¢ approaches zero where, of 
course, both F and dF dt are zero.) Another feature of the solutions of 
equation (17) is that, whereas dF /dt begins by increasing from zero faster 
than dF dt found from equation (20), it tends to its limiting value more 
slowly. Thus the approximate values of ¢ for which dF dt reaches the 
value 1-99 are 2-7 and 3-2 in the solutions of equation (17) for which R 
has the values 99-4 and 8-5 respectively, while the corresponding value 
of t in the solution of equation (20) is about 2-65. In other words the 
range of values of ¢ in the boundary layer is increased when equation (17) 
is used. While this extra thickening of the boundary layer is of more 
importance than the very small numerical changes in F(t) and dF /dt it is 
clear from the figures that neglecting the term in # in equation (17) 
produces no fundamental change in the solution, and that for even 
moderately large values of R the change is insignificant. 

In order to justify neglecting the terms involving 7 in tke left-hand 
sides of equation (12) and (13) it has been assumed that y/e > 1. If these 
terms could be ignored only where this condition is satisfied, the solution 
would not be valid at any point near the leading edge of the plate, which 
is the line » = ¢ = 0. In fact, now that the terms of equations (12) and 
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(13) can be expressed as functions of the solution of equations (17) or (20), 
it is possible to compare their magnitudes. 

The terms involving 7 are numerically greatest when 7 = 0, and in this 
case the magnitude of the ratio of the first terms of the left-hand sides 
of equations (12) and (13) to the terms involving 7 is 


20f'( f+ g) RU(1—2¢)F(d?F jdt?) F(d? F /dt?) 

rs ——TFod Fido? > (F—adFd} 
This ratio has the value 2 for ¢ = 0, and decreases to zero as ¢ increases 
to infinity. It follows that, near enough to the plate (where |), and 
therefore ¢, is small enough), the terms involving 7 will be small in com- 
parison with the first terms of equations (12) and (13) if 1—»/¢ is large 
enough, whatever the Reynolds number and however small 7 may be. 
In other words, the region in which the terms involving 7 may be ignored 
extends right up to the leading edge of the plate. The boundary layer in 
which it is permissible to ignore the terms in € and 7 in equations (11), 
(12), and (13) consists in fact of a region defined by equations like |{| < e, 
7» ¢ > A (A a suitably chosen large constant), which is a layer with a 
sharp leading edge coinciding with the leading edge of the plate. 

The boundary layer is not, however, the only region in which the terms 
in € and » in equations (11), (12), and (13) are insignificant. The second 
term on the left-hand side of equation (7) may be written as 


(arg f+g) ¢ (;2) U*atg® é fA 


th, c& 


E\h3 Sh, c&\hi, 

It is the second of these two terms which gives rise to the term in € in 
equation (11). The magnitude of the ratio of the second term of (22) to 
the first is 


g F t(d F dt) 
2(f+9) 2F 


and this decreases steadily from the value } for ¢t = 0 to zero as ¢ tends 
to infinity. Thus for sufficiently large values of ¢ (that is, at sufficiently 
great distances from the plate) the second term can be ignored in com- 
parison with the first, and in that case the term in é will not appear in 
equation (11). In fact the term giving rise to the term in € in equation 
(11) is everywhere comparatively small and for most of the plate is 
insignificant. 

Similar considerations apply to the terms in € and » in equations (12) 
and (13). 
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The cartesian components of velocity are 


ny dF U'a® {t(dF dt)— F 
“= Ll - ——_—S—n— eee 
"dt 2h t 

; oa MS 1)(1—») (t(d F dt) - F} 
~ 2h (= ( 
Ua? 


i(é 
2h i . 


> 





~¢) t 





\! (dF dt)— F; 
; 


and on the plate 


(=) | U sai] . (=) (24) 
Oz ]}..9 4 Fat\ dt}, 4 C2} 2-9 
where & — (4ax—y?) (4a) is the distance from the curved leading edge of 
the plate, measured parallel to the axis. Equations (24) imply that the 
skin friction on the plate acts parallel to the axis. The first term in the 
x-component of velocity is O(1) in the boundary layer, while the second 
term is O( R~'), the y-component is O( R~'), and the z-component is O( R~'). 
It follows that for large Reynolds numbers the flow is practically two- 
dimensional, the y-component of velocity being negligible. 

The excess of pressure over the pressure at infinite distance from the 
plate is from equation (10) 


— 1a? p/ ( , ,aF 
= Fik ] 
4h2t? | , 4 ( 


d* F\ 
dt? |’ 


2¢) 


2 25) 


Pp 
and this is O( R-') in the boundary layer. 


3. Moore’s approximations 
The solution obtained above can be compared with the solution obtained 
from the general approximate method of treating boundary layers on flat 
plates, due to Moore (1). If the plate is z = 0, 2 > d(y) and the constant 
velocity at infinite distance from the plate is (U’,0,0), the equations of 

flow are taken as 
eu 
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A solution exists with 


yy ae | ; (SVG F), 
= a i dt 


where F is the solution of equation (20), t is 4(U’/v#)!z and # x—d(y) 
is distance from the curved leading edge measured parallel to the x-axis. 
In this solution, on the plate, 


c 1’ Ri (d2F r 
(.") ~ « Fail dt? } af | (=5) 


If v(y) is taken as y*® 4a and the two solutions are compared, it is at 
onee observed that, whereas the solution derived here is fully three- 
dimensional, has no singularities except at the leading edge of the plate 
and satisfies all the boundary conditions, Moore’s solution is effectively 
two-dimensional, has singularities at every point of the cylinder y*? = 4az, 
gives no velocity components (other than (U’,0,0)) outside this cylinder 
and involves a z-component of velocity which does not satisfy the boundary 
conditions at infinite distances from the plate. On the other hand, if the 
Reynolds number is large the solutions are practically identical inside the 


boundary laver. When ¢ is small 22% = —4aZ, so that 


\(U vtyiz = RY—L(I—O}! 


and if terms of order R-! in equations (23) are ignored, the velocity com 


ponents given by these equations will be identical with those of equations 
(27). Moreover, the skin friction is the same in both solutions when the 
term involving R is omitted from equation (17). In this case, then, 
Moore’s solution predicts the form of the flow near the plate as well as 
does the more complete three-dimensional solution developed above; and 
from this fact, together with the evidence about edge effects considered 
in a previous paper (5), it seems reasonable to expect that this particularly 
simple approximate solution of the boundary-layer equations will give 
results as good as any first-order boundary-layer approximation in cases 
where the plate has a leading edge of a more complicated shape than the 
parabola considered here, and provided that the leading edge has no dis- 
continuity of slope. 

The fact that in both solutions the skin friction is the same and the 
flow in the boundary layer very similar might have been anticipated from 
Kaplun’s work (4) on the two-dimensional boundary layer. Kaplun has 
shown that in two-dimensional boundary-layer flow with external flow 
which is irrotational, the choice of the coordinate system in which the 
solution is obtained does not affect the skin friction and has little influence 
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on the flow inside the boundary layer. It seems natural to expect similar 
results with the three-dimensional boundary layer. 

Kaplun has also shown for two-dimensional boundary layers that, if 
the external flow is irrotational, it is possible to find a system of coordi- 
nates (which he calls optimal coordinates) such that, outside the boundary 
layer, the boundary-layer solution in these coordinates contains not only 
a term representing the undisturbed external flow but also a second 
irrotational term taking account of the displacement of the streamlines 
in the boundary layer. Now in the solutions of equation (17), dF /dt = 2 
and F(t) 2t+-A effectively outside the boundary layer, where A is a 
constant. Thus the velocity given by (18) is effectively 


“dl Va ws : (0.0. aA 
h, hy hy 2hat 


The velocity (1a h,, Vajhy, Va/hy) is the undisturbed external flow while 


the second term is F Uad 


[—c—g)) ) 


0,0, 
(Rh, 


and is irrotational. Thus the coordinate system defined by (1) is an 
‘optimal’ coordinate system for the three-dimensional problem discussed 
here. 


4. Conclusion 

The simple three-dimensional solution obtained here is an improvement 
on Moore's approximate solution in that it behaves better outside the 
boundary layer, and it has the advantage that the magnitude of the terms 
omitted can readily be compared with those retained. But near to the 
plate the solutions are in a very close agreement and give the same skin 
friction on the plate. It is concluded that Moore’s approximate solution 
can be expected to give results as good as any boundary-layer solution 
for plates with leading edges of arbitrary shape but without discontinuities 
of slope. 
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SUMMARY 

[n a previous paper (1) it was shown that the task of determining the displace- 
ments and stresses in a porous elastic medium could be facilitated by the introduc- 
tion of two displacement functions and by the use of repeated transforms. Two 
problems of the serni-infinite body to the surface of which a uniform pressure is 
applied along an infinite strip or over a circular area are treated here in parallel. 
Expressions for quantities of engineering significance are derived when the surface 
is either fully permeable or completely impermeable to the flow of pore water, and 


the relation to the solutions of the equivalent elastic problems is discussed. 


1. Introduction 
IN an earlier paper the authors reformulated the equations governing the 
stresses, displacements, and pore-water flow, in a loaded stratum of 
saturated clay (1). The theory of diffusion in a porous elastic medium is 
due to Biot (2); but we were able to show in the paper cited that, in 
problems of plane strain or axially symmetric strain, the dependent 
variables (the pore-water pressure, displacements, and stresses) could be 
expressed in terms of two displacement functions KE and S. The use of 
these functions in conjunction with linear transforms leads to a simple 
and practical method for solving problems relating to the semi-infinite 
body. We illustrate the technique by considering two such problems here: 

(i) The stratum z > 0 is loaded on the plane z = 0 by a uniform normal 
stress f along the strip —b < x <— b. Following Biot (2), we assume that 
the medium is isotropic, that the stress-strain relations are linear, and 
that the elastic constants and permeability are constant in space and time. 

(ii) The stratum is loaded on the surface z 0 by a uniform normal 
pressure f over the circular area p < b. 

The solutions of these problems are given in section 2 below in the form 


of infinite integrals which are in general rapidly convergent and easy to 


compute. A parallelism between plane strain and radial symmetry was 
established in (1) and this enables us to treat problems (i) and (ii) together 
with considerable economy of analysis. We have sought to simplify the 


(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 2, 1960] 
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exposition by transferring to an appendix analytical matter which is of 


special rather than general interest; in this way we hope to exemplify 
more clearly how the proposed technique might be employed in other 
problems of the semi-infinite body. 

The feature which distinguishes the problems considered here from 
corresponding problems of elastic theory is the presence of a time- 
dependent pore pressure which steadily diminishes (on the average) as 
consolidation proceeds and vanishes in the limit ¢ —- « — the steady state. 
The presence of this new variable necessitates a boundary condition on 
the pore pressure at the surface z — 0 in addition to the familiar boundary 
conditions of elastic theory. We consider two such boundary conditions 
here, namely (a) that the pore pressure vanishes at z — 0 (completely 
permeable surface), and (b) that the normal gradient of the pore pressure 
vanishes at z — 0 (completely impermeable surface). These conditions 
are, of course, idealizations of conditions realized in practice; they define 
upper and lower bounds to the pore-pressure distribution near the surface. 

The steady state distribution of stresses in a semi-infinite medium which 
is loaded at the surface was studied in detail by Love (3); the initial-state 
distribution (f — 0) is also a problem of the theory of elasticity. Both 
these limiting states are of some engineering interest, and we show in 
section 5 how some of the main results of Love's analysis can be recovered 
as limiting cases of the problems stated earlier. 

The practical aim of the present paper is the evaluation of certain 
quantities which are currently needed by soil mechanics workers in the 
analysis of field and research problems. Expressions for these quantities 
are recorded in section 7 and some results which have already been com- 
puted are exhibited graphically; more extensive computations are at 
present proceeding and will be published later. 


1.1. Historical note 

The problems stated in the introduction have been treated by earlier 
writers. Biot (loc. cit.) considered the effect of a sinusoidal load distri- 
bution acting on the surface of the stratum z >> 0 and determined the 
consolidation settlement of this surface due to a strip load by super- 
position. His analysis is closely related to the Fourier transform technique 
which we employ in this paper but is perhaps less convenient when more 
than one quantity is to be determined. More recently, Murayama and 
Akai (4) gave a solution for the pore pressure in the medium, basing their 
work on an equation of Terzaghi’s (5) which appears to be inconsistent 
with Biot’s theory. Mandel (6,7) has derived expressions for the surface 
settlement due to a point load; but, so far as is known, his work has not. 
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been extended to more general types of loading or to the determination 
of the stresses and pore pressure. The second problem mentioned in the 
introduction has been considered by de Jong (8) who used a method in- 
volving three stress functions to obtain the surface settlements. 

These references suffice to indicate a recurrent interest in consolidation 
analysis; this, despite the fact that the linear theory of consolidation 
(which is used by most workers in this field) is not fully consistent with 
the experimentally determined behaviour of clays under load. A theory 
of sufficient generality to embrace the known stress-strain properties of 
clays is an impractical quest. The practical need is to obtain theoretical 
estimates of consolidation settlements in sufficient detail to correlate 
existing settlement records, and to determine theoretically the pore-water 
pressures in loaded clay strata. Much of the earlier work on consolidation 
analysis was concerned with the estimation of consolidation settlements 
Settlement remains an important consideration in practice; but in recent 
years engineering research has emphasized the importance of determining 
the pore pressures since we have as yet little insight into the complex 


phenomena of two- and three-dimensional consolidation problems. 


2. Formal solutions for the consetidation of a clay medium 


2.1. Formal solution for the displacement functions (strip loading) 


We quote here certain results which were derived in the preceding paper 


(1), and adhere to the notations of that paper. In particular, we shall 
make use of the non-dimensional form of the governing equations (cf. 
subsection 2.5 of 1) by choosing 4 as unit of length, 6? ¢ as unit of time, 
and f as unit of stress (¢ is the coefficient of consolidation and its dimen 
sions are L?7'-'). The problem is then specified by two non-dimensional 
elastic constants: (, the ratio of the shear modulus to f, and a non-dimen 
sional ratio ». It was shown that plane strain or axially symmetric problems 
could be reduced to the determination of the functions FE and S which are 
governed by the differential equations 


" pop é2? 


where x, z are non-dimensional cartesian coordinates, p is the non-dimen 
sional radial coordinate, and ¢ is the non-dimensional time variable. The 
initial state is defined by ' 

. V2E 0. (2) 
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In the plane strain problem we introduce the successive transforms 


» \ 
E(x, z, 
rid 

i 


| E(£, 2, eos xf dé 





/ 


and 


: [ S(E.z, theos x€ dé 
7 


0 


.t) a | eMS (€,2, p) dp 
2m . } 
yY ‘f 





the partial differential equations (1) are seen to be equivalent to 


a i Ais 
"lh. = 0 
(a ‘ r)({is : )e, | 
d? : 
(i #)s, . 


(5) 


Choosing solutions of (5) which vanish as z approaches infinity, we may 


E. Aye of 4 Aye af?+p)! 


€ 


S. = Bye ’ 


write 


(6) 


where A,, A,, B, are functions of p and € which are determinable from 
the boundary conditions. 


2.1. (i) The boundary conditions of the formal solution (strip loading) 


We shall first consider the case when the surface z = 0 is impermeable 
to the flow of the pore water from the underlying stratum. The boundary 
conditions of this problem can be stated as 


Ca - 
-=O=— onz=O0fort>0,—-w<r<om, (7) 


Ce 


jt. jal < 


and o.. = 
10, ja! 


zs 


onz = 0,t > 0, (8) 


where o,,, ¢,, denote the shearing stress and the compressive stress, and 
a is the excess pore-water pressure; (8) prescribes a uniform unit load on 
the strip 2 <1. 





214 J. MCNAMEE AND R. E. GIBSON 


It was shown in sections 3 and 4 of (1) that o,,, o are defined in terms 
of £, and S, by the repeated integrals 


7? eo Be 1B. _ dS, 
| Esin(xr€) dé —, | on — +2 | dp 
7 ~7l az az 


. 
0 Br 


a 2 7 ] 5 d? ad ds. 
| cos(ré) d& - ert 7 _--é le ‘ dp 
2G. C1 2m . dz* dz 


Ded 


=Ti 





0 Br 
In these expressions ( is the shear modulus when the unit stress is f and 
7 is an elastic ratio; the symbol Br is used to signify the Bromwich- 
Wagner contour y—ix to y+im%. 


If o,, and ¢o @z are to be zero on z = 0 for all x and for ¢ > 0, we must 


have dB 
0 
dz | 
Z 0; 


\n @S 
and £2) E : o| 


dz 


i.e. using (6) ¢ 


fA, +-(€24 pyAa, —s 0) 
np(e?+p)tA,+é*B, 0 | 
The compressive stress o., is defined by (1, sections 2.4, 3, and 4) 


» . | + 
cos(.r&) dé — | 
T ' 2a 


. . 
0 Br 


.. 


2G 


ore +(4 208] ap 


and the boundary value of o,, can be written as (cf. (6) and (10)) 


© 
. 


a it ¥£ . aoe 
22 n<* 27 23 2 2.) > 
Eat ge | el F [£3 -+-(np—&*)(€2+- p)* |eos x€ dédp. 
; 0 br ‘ 
The boundary condition (8) can be satisfied by means of the known dis- 
continuous integral 
Sin € COS(XE) | oy dédp = jt, (11) 
pe \0, |z| > 
which is valid for ¢ > 0. The constants A,, A,, B,, are now given by 


(l+s)t . \ 
Fal sin € 


0 Br 
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if we introduce the change of scale 


p= fe (13) 


in the p-plane and use the abbreviation 


A = 14(ns—1)(1+8)!. (14) 


2.1. (ii) Completion of the formal solution for E, and S, 

Expressions for the pore pressure, stresses, and displacements, are of 
little interest at this stage unless we wish to evaluate one of them inde- 
pendently of the others. We proceed to determine the functions FE, and S, 
defined in equations (3) and (4). Using the results of the preceding section, 

sin € 


2GE, = [ - tei as, (15) 


2mé J sA 


Bs 
sin€ [ e&€*s-€(1+-s)! . 
—% Dmié? —- en _ ds, (16) 
” Bs 


Bs denoting the Bromwich-Wagner contour in the s-plane. 
The denominator A can be written in the form 
ja ale 2n—D/o 
: (ns— 1)(s+1)!—1 
aa (re t40 22) 


2n 











where 


2) oa i 

7 7 
The zeros of \ are s = 0 and s = —A where J is always positive and less 
than unity since 7 is finite and greater than or equal to unity. 

The contour integrals for E, and S, reduce to circuits around the zeros 
of A, together with a loop integral (extending to —oo) around the branch 
point s = —1l. Noting that E, has a double pole at s = 0, we easily find 
that 


and 


sin 


2GE, =~ Jere vy, —14Y v4 


2n—1- 
2e—i—w)yr 
+ Toe peel 
L— Qwte—l-er 
(n+1)?—(2—nw)* 





“wt —ate-s]| 
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where + = t&?, y = 2, w = 1—A, and 


] 2 [ p _{ 1+ n(v? + 1) |v? cos vy— v sin vy 
=-—; —— 


(v?-+ 1)*(v? + w){v? + (1/n2w)} 





eT? dv _ Ee | 
(v? + 1)?(v?-+-w){r?+(1 n*w)} aa dy —_ 


€ Tr 4,2 dv 


» bd 
ait) FE Fay T OD 





The expressions for these integrals are only of interest for reference and 
they are recorded in an appendix. The stresses and displacements are 
expressed as derivatives of EZ, and S.; but it is not very convenient to 
obtain them by forming derivatives of the expressions (18a) recorded in 
the appendix. It is usually much simpler to evaluate directly the deriva- 


tives of (18), e.g. by using the known integral 


e-2" cos yv * «) 2» y es a 2 
— dv e%8*| e-8y erfel Bai — 7 |+ Fy erfe| xy) + — 
vy? +B? 4p | 2at Ba! + 2a 


. 


(B + 0, a > 0). 
The final solution for the functions E, and S. is now complete. So far, 
the analysis applies only to strip loadings, but it will be shown in the 
next section that the extension to radially symmetric loading is trivial. 
2.2. Formal solutions for the displacement functions (circular loading) 
It was shown in (1) that in the case of radially symmetric strain E and 
S are unchanged if we interpret the operator V? as 
vr <a l oe F< a 
cp” pcp cz 
p being the radial coordinate. Instead of (3) and (4) we use the transforms 


E(p,z,t) = [ E,(€, 2, &dg(p&) dé 
0 (3’) 
E,(€,2.t) = a | By é.2.p) ap | 


7t 
. 


Br 


and . 
S(p,z.t) = | S,(E, 2, HETy(p&) dé | 


+ 


eS, (€,2, p) dp | 
2m . 


Br 
the suffix / indicating a Hankel transform (see 1, section 3). 
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The functions EZ, and S, satisfy the same formal equations as £Z, and S., 
i.e. equations (5). Apart from the replacement of Fourier by Hankel trans- 
forms, the only new step is the satisfying of the boundary values of ¢,, 
on the surface z = 0. Since o,, is required to have unit value inside the 
unit circle and vanishes elsewhere on z = 0, we use the discontinuous 
integral 

- as “ 
sa | | Nee Wh(6)™ dpdg = ‘ eS) UDO. 
0 Br 
in place of (11). Equations (6) still hold when the suffix c is replaced by 
h and the constants A,, A,, B, are given by equations (12) if sin€ is 
replaced by J,(é). The expressions for E, and S, can be written down 
from (17) by writing J,(&) in place of sin €. 


3. Verification of the solution 

The solutions of the problems considered are defined by the pair of 
functions £, and S, or by the pair E, and S,. It is desirable to verify, at 
least in outline, that these functions behave as we should expect them 
to do on physical grounds. For example, it is known that the normal 
displacement w is infinite in the strip-loading case. It is simple to establish 
directly the convergence (or divergence) of E,,, and S., and their relevant 
derivatives, as we show briefly below. 


3.1. Convergence and differentiability of the integral representations of E 
and S 
It will simplify the discussion to introduce a notation which covers the 
plane strain and the axially symmetric solutions. We shall write 


(19) 


K(. é) " 7, cos(x€)sin é 


’ Jol pe (€) 


and adopt the convention that the upper line is to be used in the plane 
strain solution and the lower in the axially symmetric case. The ex- 
pressions for E and S can now be written as 


r 1 x 1+-2 
2GB(*,2,t - ta (F d |- TY ¢-y 
p ) 2e Ya 2-1 * 
2e——w)r 


+ (—el[@ +l —2—7a)'] 





(ewe — aoe») eth eM), (20) 
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Zz 2 
é) v dé ei, 
p 2n—1 
2w! , 
e~' 
(7 + 1)? —(2— nw)? 


wr). (21) 
It is clear from equations (I8a) that convergence at the upper limit is 
easily verified when z, ¢ are positive or zero. 

3.1. (i) To investigate the behaviour of E at the lower limit we expand 
the expression enclosed in square brackets in equation (20) in powers of €, 
and retaining only the term independent of € we obtain 

2(1—wt) 


> 


(1—w)| (7 + LD? —(2—yw)?] | 


yw  Pu(l— nw re) 


w)? (1—7*%w)?* 

w yw? I 1) Jert 2 Oe) 
1—*w*\(1—7*w)® (1—w)* 2tt 
or possibly O(€), and it can be shown that the above expression is 
identically zero. From (20) and (18) it can be established that the limit 
is, in fact, O(€*). 

Since A(x, €) is O(1) and K(p, €) is O(€), € — 0, it is clear that the integral 
representation of E converges at the lower limit for finite ¢. An additive 
constant is needed when ¢ -- «% and we shall consider this point in section 5. 

The convergence of the expression for S is determined by the behaviour 
of the integrand at the lower limit. For €—- 0 it can be shown that the 
bracketed expression in equation (21) is O(1), and consequently the Hankel 
integral for S is convergent but the Fourier integral diverges. However, 
by adding a suitable constant, e.g. 

2sin€ 
ne? 
to the Fourier integrand convergence can be secured for finite z and t when 
7» = 1. The anomalous behaviour of S in the limit ¢-> © is well known 
in the corresponding elastic problem. 

We can easily confirm by inspection that derivatives of E and S of all 
orders exist when z is greater than zero. In determining the existence of 
the boundary displacements, we observe that for ¢ > 0 we need consider 
only those terms in equations (20) and (21) which do not possess the 
dominating factor exp(—r). It is then clear that 

COE cK ek @E @E BE 
cp” pcp” ox’ ep?” éx2’ Gz? 
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all exist on z = 0 since the infinite integrals 


J p 


| K(*.¢) dé ond | K,(*.¢) dé 
’ iaeens é 
exist, A being defined by (19) and A, by 
= sin(é)sin ¢ 


k, ( J "y 
‘ hlpe rhe) 


Similarly @S éz exists on z = 0. Finally, the integrand of 
ek 
Op C2 
vanishes identically when z approaches zero. These considerations ensure 


that the solution can be completely verified. 


4. Solution when the boundary : — 0 is fully permeable 
The analysis here is very similar to that of the preceding case and we 
record only the final expressions for E and S. The boundary condition 
on the pore-water pressure is now 
o=0 on z=090 
and the considerations of section 2.1 (i) apply mutatis mutandis as before, 


We tind 


, r 1 oe .3 

20k(*.2.4 K( P T(I,4+-e “I,)—e-¥ 

p | - Pp ji = : 
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208(".2.4 ; nf -&(8)| a “1 v dé, 
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0 
where 


(e+ 10? + (y—1/n)*] 


oo v? cos(yv) + | n(v?+-1) I |v sin(yv) a 
a 
(v?-+ 1) v? +-(»— /n)?] 


, v3 d 
eerie. »—TY v 
= (v?-+-1)[0* + (4 —1/n)?] } 


The integrals (24) are similar to the earlier integrals (18), and they are 
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given in terms of computed functions in the appendix equations (24a). 
The forms (24) are preferable to (24a) in analytical derivations. 

The convergence of the expression for E at both limits for finite values 
of t can be established by the argument of section 3.1. However, as before, 
to secure convergence of the Fourier integral for S requires the addition 
of a constant to the integrand; e.g. 

2 sin & 
ane? P 
and S becomes 


6 p » 
2GS(zx,2z,t) = = 2 -ms = ie "Ig}e-¥ cos.xé dé 


which is convergent for finite values of z and ¢. 


5. Limiting forms for the functions E and S 

It may be useful at this stage to consider the relation between our 
solutions of the diffusion problems and the solution of the equivalent 
elastic problems. The steady state solution of the diffusion equations we 
have considered is identical with the solution of the corresponding elastic 
equations since the excess pore water pressure reduces to zero, as t > &. 
It is evident from equations (20) to (23) that as t > o the stress functions 


become 


K(*, ee +y) dé, 
¢ \P 


2 flelx 
IGS. = - l ~ KI" ,é\e-" dé. 
z—1 | eK (hee . 
0 


The convergence of E can again be secured by the addition of a suitable 
constant; the displacements and stresses are not affected by this addition 
since they depend on space derivatives of E. The Fourier integral for S 
diverges at the lower limit but we can, as previously noted, secure its 
convergence for all finite values of z. The vertical displacement w is given 
by (1, section 2.2) 

w= 


and it is evident that the addition of a constant to S is equivalent to 
a rigid-body movement of the whole medium—the horizontal displace- 
ment and the stresses are obviously unaffected by this movement. We 
are then unable to ensure that w in the plane strain problem shall be 
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everywhere finite for all values of ¢ but it is perhaps permissible to ascribe 
some physical significance to the difference 


wW—Uy 9 


which is finite for finite values of t. This difference is usually termed ‘the 
consolidation settlement’. 

We can obtain the elastic displacements and stresses directly from E, 
and S, by differentiation, but we shall confine ourselves to the derivation 
of the surface settlement w,_, and the dilatation ¢ in the axially symmetric 
case. The surfaice settlement is given by 


nf APES) ge 
(Qy—-NG)~ & 
0 


This integral can be expressed in terms of complete elliptic integrals of 
the first and second kinds of modulus p or 1/p (see, e.g., Watson (9)), 


2n 7 
= 0 < =< 
2 1) Ele) ( p ) 


5) 
—. |(; (1 =)&(-)] (p > 1). 
7G(2y—1)| \p p*} \p, 
which was given by Love (3).+ The final dilatation is given by 


VE Ip(pé (Ee dé. 


l [ 
: Gi(2y—1), 
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A somewhat different expression for e in terms of incomplete elliptic 
integrals (modulus /) can be obtained from Love's analysis (loc. cit.) in 
the form 


pyle aK Bu) HK E')—K'ksnp] (0 <p <1) 
7) 7" 


Qian = 


l 
7G(2 
G2n— pe Sele T- K’ksnp| (p > 1). 

Tr =7) 
where K’ and E’ are complete elliptic integrals of the first and second 
kinds with modulus k’, and where 
z = kr,sn(p, k), 1+p = r,dn(p, hk), 
1—p| = kr, en(p, k), rz; = 27+(1+-p). 
This suggests that the Bessel function integral is in fact equivalent to 


+ This slight and regrettable ambiguity of notation was not suspected at the outset. 
The elliptic integrals are not encountered elsewhere in the paper, and modification of the 
displacement function and kernel notation is scarcely justified. 
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Love’s expression, and we can establish this equivalence by specializing 
a result of Watson's (9, 390). 

Similar correspondences with known solutions may also be established 
in the case of the strip load. It may also be noted that initially (t = 0) 


r 


‘ é}; v dé, 
p 


, : l * ; 
E--0 and S + ys) | Al 


which are identical with the values of E and S in the limit as 7 - oc. Our 
solution for t = 0 (before any dilatation has had the opportunity to 
occur) is therefore seen to be the elastic solution with the restriction of 
incompressibilit y of the medium (7 = «); this result could, of course, have 


been anticipated by physical reasoning. 


6. Displacement functions for the point load 

\ further limiting case of some importance is that of the point load. 
The displacement functions for this case can be obtained from equations 
(20) to (23) by allowing the radius 6 of the circular pressed area to approach 
zero while the total load P = zb?f remains constant, but in order to carry 
out this process we are forced, temporarily, to abandon dimensionless 
variables. We introduce the dimensional quantities 

E’'=E, S’=b6S, ¢ =O%/c, £ = €/b, G’ = fG, ete., 

into equations (20) to (23) and then drop the primes. With f = P/zb? 
and 4 —» 0 we find that the formulae for EZ and S appropriate to the point 
load are given by (20)-(23) with A(p,&) replaced by 


EJ, pé ). 


La 
=7i 


From these displacement functions formal solutions to a variety of prob- 
lems of the semi-infinite porous medium, under vertical load distributed 
in some manner over the plane z = 0, can easily be obtained by super- 
position. 

An alternative and perhaps simpler treatment of the point load case is 
to make use of the repeated Fourier integral solution in (1, section 7). 
The equivalence of this with the Bessel function solution is easily estab- 
lished using Poisson type integrals for Bessel functions. 


7. Expressions for quantities of engineering interest 

The quantities of most interest to engineers are the surface settlement 
W,_9, the excess pore-water pressure co, and the dilatation e. The predic- 
tion of the time settlement of foundations on clay strata has hitherto 
been almost entirely restricted to cases where the diffusion of pore water 
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is one-dimensional. However, problems frequently arise in practice where 
the assumption of one-dimensional flow is not even approximately satisfied 
and where an assessment based upon a one-dimensional theory could lead 
to a significant underestimate of the rate of settlement of the foundation. 

Variations of volume and of excess pore-water pressure in the clay lead 
to changes in its strength, and some knowledge of the pattern and progress 
of these alterations in the vicinity of the loaded area can be used as a guide 
in assessing whether the foundation can safely carry additional load. In 
this section we derive explicit expressions for these quantities. 

The displacement w of the medium is given (1, section 2.2) by 


cE 
On the surface z = 0 the first two terms vanish and we have 


4 < 
Hi; aks (?. é) FigM) ae. 


where F denotes the integrand appropriate to the impermeable or per- 
meable boundary z = 0 (equations (21) and (23)). 

In both cases this integral can be shown to possess a discontinuity at 
t = 0 corresponding physically to the immediate settlement of the surface 
of the medium upon application of the load. This immediate settlement is 


Ws-0t=0 = Mi fi é Kt’ é) dé. 


For the Bessel kernel this has been evaluated already (section 5). To 
the Fourier integrand we may add the constant —sin£ 7G&é* to secure 
convergence; thus 


[ w(0)— w(x) ], 94-9 = 1+2 log 1+.2x + 1—r log 1—2x)], 


_ 
G 
which is a known result. 

The additional settlement (consolidation settlement) will be given there- 


fore by 


‘x 


a 7; L it 
2G[w—wy, o|.-6 = ee [ -*(7-¢)]ertier— 


. 
) 


€ ~{(27- 1/77" erfe ] l ee] dé (25) 
7) 


when the surface of the medium is permeable. 
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This integral may be expressed in terms of tabulated functions in the 
special case 7 = | (Poisson’s ratio of the medium zero) for the strip load, 


. t\} ltox fl—2z 
2G[w uy ol: 0 (")° erf (= r} L_erf ( a —_ 


I wash (1-+2)*) 
lt+a)E -}+(1 
al *) ‘| 4¢ | ( 


and for the centre of the circular load 


(1 t\)} 
2G[w—uy,ole-0 erfe( 5) } 2(< 


For the impermeable surface 


x 
7 


— | 7K ("2 erf(ét#)+ 


=) a 
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2G wu el . 


v(2y—1 
a 1 = ) he 0-08 1 Lert (Eeit#)}] 
+ Uf ‘ yw 


2n—1 set th 

_—— (; ,-1-anwntterfe(S “\\| ag. (26) 
m* n(n 1)w'| nw | 

These integrals converge in all cases for finite time, but as ¢ approaches 


infinity the Fourier integrals diverge. 
The dilatation ¢ is related to the function EF alone by 


f V7E, 


and we find 
] — -- in 2 
20 | K(".¢) f ‘erfel étt) 4 
2n—1, — 
i 


7 l {(2m—1) 92 }E*t+(y inks erfe( eps _ 
») ” : 


for the permeable surface, and 


] aE 2 
=> | K(*.¢] ¢ € erfe| &tt) — 
27 l , p 20 
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9 i » 
(27 —1)w* p~(1—w)€%t fe erfo 5 fait!) 


D4 hl 
HY? 
~~") d& (28) 
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for the impermeable surface. 
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The excess pore water is given in terms of the quantities e and @S/éz 


(1, section 2.2) by the equation 


os 
oO 24/ ne}, 
Cz 
where 


*.¢). | 1--erf(ét!) 4 


] l etn varieterfe™ , “é0) | dé (29) 
” ”) 


for the permeable surface, and 


, 5) 
; .é) | ; erf(ét!) e-l ott | erf (Ew! t*) | | 


2n l 


etl anroneterte(§ : — 
mn 3—n(9 + ew} a” 


for the impermeable surface. 

The consolidation settlements of the centre of the circle and the strip 
loads have been computed for a medium with a Poisson's ratio of zero 
(» 1): the results are shown graphically in Figs. | and 2. More extensive 


computations are proceeding the results of which will be published later. 


APPENDIX 


The integrals /,, /,, and /, of section 2.1 (i) can be expressed in closed form as 
follows: 
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SUMMARY 
This paper discusses the problem of the deformation of a thick, isotropic, elastic 
plate containing a spherical cavity. The plate is considered to be infinitely extended 
in two dimensions and the centre of the cavity to be midway between the faces. 
It is supposed that the surface tractions specified over the boundaries are continuous, 
axially symmetric, and torsion-free. The resulting strains are assured to be every 
where sufficiently small for the classical theory of elasticity to be employed. The 
solution is found by a method of suecessive approximation, illustrated by the 
following examples 
i) Uniform tension applied at infinity parallel to the faces of the plate, the finite 
bounding surfaces being stress-free. 
(ii) A flexural couple applied at infinity, the finite bounding surfaces being stress- 
free. 
Comparison is made with the elementary, approximate solutions obtained by 
neglecting the effect of the faces. 


1. Introduction 

Tue problem of a thick plate of isotropic, elastic material, extending to 
infinity in two dimensions and subjected to arbitrary boundary tractions 
and certain internal loadings, has been dealt with extensively by Dougall 
(1) who also considered by a similar method the infinite circular cylinder 
(2). More recently, Green (3) has developed the analysis for stress systems 


in a thick plate containing a circular cylindrical hole; this was done partly 


in order to check the assumptions of the theory of generalized plane stress. 
Later (4) he solved the problem of a uniform distribution of force along 
a line perpendicular to the faces of the plate and acting parallel to them; 
also, (5), the general problem of the equilibrium of plates and cylinders 
in which the faces or plane ends are stress-free. Using Hankel transforms, 
Sneddon (6) has examined the stresses produced by pressure on the faces 
of a thick plate. The effect of a spherical cavity in a circular cylinder has 
been examined by Chih-Bing Ling in the cases of torsion (7) and simple 
tension (8). In this paper a similar investigation is carried out for the thick 
plate.t 
+ A referee has also drawn the writer's attention to (18). 


(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 2, 1960) 
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For simplicity, the centre of the cavity is taken to be midway between 
the faces and the only boundary conditions considered are axially sym- 
metric. 

In problems with axial symmetry and uniform temperature distribution, 
it is convenient to distinguish between two different types of stress system. 
Firstly, that which arises from a displacement field which, at each point, 
has direction in the axial plane through the point, known as the torsion- 
free system. Secondly, that which arises from a displacement which, at 
each point, is perpendicular to the axial plane, known as the torsion system. 
A full discussion of these systems, together with the thermal stress system 
arising from a non-uniform steady temperature distribution, will be found 
in (9).+ This paper describes the analysis for problems involving only 
torsion-free stress systems. The discussion of systems with torsion and 
thermal stress will appear in a later paper. 

The inherent complexity of thick plate problems is increased by the 
addition of an internal boundary of a shape not geometrically compatible 
with the faces. Although methods are available for satisfying boundary 
conditions on the faces and on the spherical boundary separately, it does 
not seem possible to satisfy the conditions on the two boundaries simul- 
taneously without the use of infinite series, the coefficients of which are 
determined by a system of linear equations involving an infinite matrix. 
Because of this, the iterative procedure described below seems to be simpler 
for application to specific examples. 

Similar methods of successive approximation were used by Howland (10), 
and Howland and Stevenson (11) in dealing with problems in the analogous 
two-dimensional region, the infinite strip with a circular hole. 


2. Method of solution 
Use will be made of spherical polar coordinates (r,6,y) and cylindrical 
polar coordinates (a, y,z), these being related as follows: 
rsin@cosy = weosy, rsin@ésiny = wsiny, rcos@ = z, 


The faces of the plate will be taken as z = -+-A and the spherical cavity as 
r= pp. 

When the body force is zero, the equation of equilibrium, in terms of the 
displacement vector D, is 


VivV.D) T (1 wa 2y)V°D — 0, . 1) 


where 7 is Poisson’s ratio. The corresponding stress vector R,, across any 


+ A thesis by the writer for which the Ph.D. degree of the University of London was 
awarded. 
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element of surface in the material with unit normal vector n is given by 


n j » 


R “en (V.D)n+2u(n.V)D+pn (VD), 
” 


where yp is the rigidity. 

Thus the problem consists in finding a solution D of (2.1) such that the 
associated stress vector (2.2) across each element of the bounding surfaces 
2 th, r = p, takes on prescribed values. 

A series will be constructed in the form 


D,+D,+D,+-D,+..., (2.3) 
where the terms D, (r = 0,1,2,...) are solutions of (2.1) satisfying the 
following conditions. 

D, is such that the corresponding stresses as w tends to infinity take on 

the values prescribed in the problem. 

D, is such that the stress vector R, corresponding to D,+D, takes on 

the prescribed values over r = p. 

D, is such that the stress vector R, corresponding to D,+-D, +D, takes 

on the prescribed value over z th. 

D,,., (r = 1,2,3,...) is such that the stress vector R, corresponding to 
D,,+D 
op.g (Y = 1,2,3,...) is such that the stress vector R, corresponding to 
D,,.,+ D,,.,. vanishes on z th. 


r 


,; Vanishes on r = p. 


2r 


D 


The first »+1 terms (n > 0) of the series (2.3) will provide the solution 
to a problem in which the conditions on one boundary are satisfied but not 
those on the other. The error occurs on the spherical boundary or the 
faces according as n is even or odd. If this error tends to zero as n > &, 
and provided that the series, together with its first and second partial 
derivatives, converges, and the terms D, (r > 0) do not affect the conditions 
at infinity, then the sum to infinity of this series is the required solution of 
the problem. 

In practice, the convergence of (2.3) is not required and will not be 
considered: it is sufficient that the error on the boundaries should tend to 
zero. The series is then terminated after a finite number of terms as soon 
as the desired degree of accuracy is attained. 

The torsion-free system may be dealt with by using the method of 
Boussinesg (12) and introducing two harmonic stress functions, or by using 
the single biharmonic function of Galerkin (13) or, again, by the general 
biharmonic vector of Papkovitch (14). However, the advantage of using 


the single stress function approach seems to be outweighed by the fact that 


it is biharmonic and leads to more complicated expressions for the stresses. 
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The harmonic stress functions approach will therefore be used here. In 
determining D,,,, use is made of stress functions with singularities only at 
the origin, and for D,, use is made of those stress functions with singularities 
only at z = o. In this way the conditions on the plate for large w are 
unaffected and, in numerical examples, it is found that disturbances caused 
by prescribed stresses on the boundaries are soon smoothed out and the 
error after (n+-1) terms of the series (2.3) tends fairly rapidly to zero as n 
increases. 

It will now be shown how to construct each term D, (r > 0) from its 
predecessor D, , and the value of D,, which is assumed known. It will 
also be assumed that D, does not contribute to the stress on the faces. 
For simple states of stress at infinity, there is no difficulty in finding such 
a solution D, as the later examples show. In any case, this restriction is 
introduced merely for simplicity. Even if Dy produced stresses on z = +h, 
these could be annulled by a modification of D,. 


3. Construction of D,,,, 

As shown by Boussinesq, and used extensively by other writers, the 
solution of (2.1), in the case of an axially symmetric, torsion-free stress 
system, may be represented by the sum of two solutions of the types 


2uD = 4(1—n)dk—V(z¢). (3.1) 
2uD — Vu, (3.2) 


where k is a unit vector in the z-direction and ¢, % are harmonic functions 
independent of the azimuthal angle y. 

Let @,, and, denote the values of D given by (3.1) and (3.2) respectively, 
when both ¢ and & are replaced by r”P,(cos@) (n = 0, +1,+ 2,...) where 
P. (cos @) is the Legendre polynomial of degree n, i.e. 


Zu, = 4(1— nr" P,n—V(er"P,), (3.3) 
up, = Vir"P,). (3.4) 


Now suppose that D,, may be expanded in a series, convergent within 
some sphere centre the origin and of radius p, > p, in the form 


D,, = z h meh ep, mip, +1}, (3.5) 
n=l 


where /!?"), m'?” are constants, which, for convenience, have been made 
dimensionless (apart from a factor «) by the introduction of the factor 
h n+1 


Firstly, the stress components #7, r@ of R, corresponding to D,, will be 
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found. Then a method will be given for expressing the term D,,,, in the 
form : 
D,,.1 Zz Aeth{ g(t + Dep _ | qirt Mp n+1} (3.6) 
1 


n 

such that the displacement D,,+D,,,, produces zero stress across r = p. 
In the construction of D,, a further term must be added in to account for 
the prescribed stress on r = p. However, if the prescribed stress com- 
ponents are expressed in series of the form (3.7), (3.8), the additional term 
may be found by the same procedure as given below. 

After some algebraic manipulation and use of the recurrence relations 
for Legendre polynomials, the stress components #, r corresponding to 
D,, are found to be 


<, | 2n+3 


n 


= (n-+-1)(n? —n—2—2n) al(@r) 
S| M2? 4 


n(n —1)(n—4-+-4n) 


2n—1 


An-21'2"), —n(m —1)A" 2m (2r),| p (3.7) 


jen 


x 


* {(n?-+2n—142n) yynan 
> A swig T 


a 2n+-3 
n 0 


(n— 1)(n—4-4-4n) ‘ini —" bs 
os iT" t dn-21(27), —(n—1)A" *mitr,| Pisin 8, (3.8) 

where A — rh and the dash on P, denotes differentiation with regard to cos @. 

Sternberg et al. (15) have solved the problem of an axially symmetric 
stress distribution over the surface of a spherical cavity in an infinite 
medium. Their method may be employed to find D,,,, in the form (3.6) 
such that the stresses (3.7), (3.8) are annulled. 

After some algebra it is found that 

ars 1 \2n + 1p(2r) ” \2n—17(2r) mY2n— Ly) (2 
Pr Fe Pn Ao’ l, ; 1 Pn Ao" lr 2 it Pn o my "1 


2r+1 ‘ y2n + 1p(2r) ” \2n —1)(2r) \2n —3)(2r) 
Un Un ro I +1 In ro I ee GnXo Ab 37 


iv }2n —1,,,(2r) v )\2n —3,,,(2r) 
tan AG Mri +nAo Mn—3 


as 


where A, = p/h and 
Pp, = (2n+1)(n+-1)\(n—1)e,, Pp, = (2n+1)(n— 1)(n—4+4-4y)e,, 
p, = (2n-+-1)(n—1)(2n—1)e,, dn = (n—1)(n+ 1)(n+3—4y)e, 
(2n +-1)(n—1)(n—4+ 4n)(n+-3—4m)e,,(2n—1)-!4 
+-{2n°— Lin* + 40n3 — 45n? + (26—8n?)n— 12(1+- y)(1—7)} X 
XK €,2n-(2n—1)-! 


(n—2)(n—3)(n—6+-4y)e,,_2, q's i 


n—1)(2n+-1)(n+3—4n)e,, 





(n—2)(n—3)(2n —5)e 


( 
e ¢. = hin?+-(1 2n)n-+1—n} 1 


n - n 


(3.10) 
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4. Construction of D,,,, 

Sneddon ((16) 468) has given a method for determining the stress distri- 
bution in a thick plate subjected to axially symmetric boundary tractions 
using Hankel transforms to derive the required biharmonic stress function. 
Here the corresponding analysis will be given for the two-function ap- 
proach of Boussinesq. 

The stress components 22, Z@ corresponding to the stress function 


$ = | (A cosh z+ Bsinh £)Jy(€er) dé (4.1) 
{ €%2(1— »)(A sinh €z-+- B cosh £2) - 


— &z(A cosh z+ B sinh &z)!.J,(Ear) dé 


Zs 


( £7|(1 —2»)(A cosh €z + Bsinh €2) 





&2(A sinh €z + Beosh €2)}.J, (Ea) dé 
Those corresponding to the stress function 


( &(C cosh 2 + D sinh &2z)J,(Ear) dé 


0 


| &(C' cosh &2z + Dsinh &z)J, (Ea) dé | 


(4.4) 


{ €(C'sinh £2 + D cosh €2)J,(Ear) dé 


0 


If A, B, C, D are determined so as to annul the stresses 22, Za produced 
by D,,., on z th, the stress functions (4.1), (4.3) may be taken to define 
D In the case of D,, a further term must be added to account for the 
prescribed stresses on the faces. However, if the Hankel transforms of 


ar+2: 


these stresses are found, the additional term is easily found by the procedure 
given below. 
It is convenient to separate D,,,, into two parts as follows: 


n 


zx 
~ B2n+1f y(2r +1) | gf2r 4d) ea 
D,,.; 2/ {Pon Pan +9an penal 


7 
~ fens y,(2r+ 1) (2r +1) 1 & 
t 24 "Pen 1 i) Qne17 Yon 1 y 2n+2s° (4.5) 
n 


The first part gives displacements which, in the region z > 0, are the mirror 
images in z = 0 of those in the region z < 0: this part will be referred to 
as the symmetric part. The second part, with the opposite property, will 
be referred to as the antisymmetric part. 
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Using (2.2), it may be shown that the stresses 
sponding to the symmetric part of D,,,,, are 
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0) should vanish over z 
form (4.2), (4.4) with 


n)eosh Eh 


+-h, the stress 
and B, C satisfying 
Eh sinh Eh} 


functions defining ‘he symmetric part of D,,,, are to be expressed in the 
A=D 
BE*}2(1 


C'& cosh £h S *n gne-th 
: i kant 


C& sinh fh = ’ 


> B,, EN" —£fh 
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(3.5), the following expansions are used: 


27)sinh £h+-€h cosh Eh} 
In order to express the values of 4, & given by (4.1), (4.3) in the form 


(4.10) 
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Hence the stress functions defining the symmetric parts of D,,,. may 
be written r 
db 2,4 Safle + BypSnsel Pe as 
7 (4.13) 
yp >" 2nyy 2+ Dp2n+2pP 
n 
where 
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Solving (4.10) for B, C, and substituting from (4.7) for «,, B,, (4.14) 
reduce to 
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(4.16) 
The integrals /,, J, have been tabulated by Howland (10), and their 


asymptotic values shown to be 


I s! J s! 
«~ on’ va ~ Saath 


The antisymmetric part of D,,,, may be dealt with in a manner which, 
apart from the occurrence of divergent terms in the integrals, is similar to 
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the above. The stresses 22, 
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hin this case may be shown to be 
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view to constructing the antisymmetric part of D,,,, (in which 
0) to annul these stresses, consider the unintegrated forms of 
s funetion 
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,, (4.19) may be solved to give 


(2 


2r+) 
1 


2)! 


,,(2r 
2th 2m ~Yam 


(2m 


| 2£h] (4.20) 


(2m 


(2r +1) 


2m 1 


| —2y)— 2een2] 2” 
M 1) r ’ 2)! 


2m 
qi2r 1) 
| —(3—4m)-+ 2h-+<¢ fh] lami | 


‘(2m —-3)'f? 


(4.21) 


f cosh &2J,(€a) and &sinh &2J,(€o7) respec 


tively, the expressions for 


A and PD could be integrated with respect to 


from 0 to », these integrals would then give the required stress functions 
¢, & defining the antisymmetric part of D,,,». Unfortunately, the values 
of A and PD) given in (4.20), (4.21) contain terms which cause the integrals 
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to diverge at the lower limit. Expanding the integrands (4.18) in powers 
of €, these terms are found to be 


(2r+1) (2r +1) 
a(1—n) Pt 4s 31 — ny 


2p 4.22 
he * & oa he a 


yr, h rP, 
, é 


| PP, (4.23 
hé é mp J 


Po wl 
(1 —»)(1—2») 2h 


3(1—)(1—2y)py"’ 
respectively, where 
f= 34 2y)p PY + 2 ype) — dag”, 
g = fl 2(1— n)(1 2m) +5] pyr? — BCL 9) 2) pg + AL 2 grt. 


The first two terms in (4.22) are independent of z, w and the first two 
terms of (4.23) contain only the factor rP, (= z). These terms, therefore, 
give rise only to a rigid body translation and so may be removed from 
the expressions (4.18) without affecting the stresses. The terms of order 
& ' will, however, first be multiplied by the factor ¢ ‘ in order to secure 
convergence at the upper limit on integration. 

It is easily verified that the last term in each of (4.22) and (4.23), taken 
together, do not contribute to the stresses across 2 th. These terms 
also may be multiplied by e * and subtracted from their respective ex- 
pressions (4.18) which now become integrable. 

Each term of (4.22), (4.23) is a harmonic function and so the required 
stress functions defining the antisymmetric part of D,,,, are 


d, | eA cosh &zJ,(Ear) 


. 
0 


+1) 


yp | |ED sinh 2y(Ee) + 3(1—7)(1 2) hes 
| (2r4 Vy th 


+ (1 2y)(1 7)? 2h (23 jewt)} dé (4,25) 


since r?P, == 22—jo?, PP, = 23— Fz". 

The significance of the divergent terms in the integrands will be investi- 
gated later. 

Proceeding as before, and omitting terms which give rise only to a rigid 
body translation, the above stress functions may be expanded as 


b= Shemp, a Shige, (4.26) 
n nel 
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where 

— % sinjz (2r+1) ) Qn (2r+1)) 

= > t " Eom 1P2am-1 1 a EY | -S8 
m=1 


zx 
» faney (2r+1) 4 2 (2r+1)) 
> \ "Gom—1 Pom —1 "Hom, 192m—15 | 
m=1 
in which 
l 
one (ye , 197" ¥ ’ 
"E, 2)'(2n)! (3—4n) lon .2n-2 2Tomsen—1 t J oms2n 2s 


» 


3)'(2n)! 2m+2n—-2 


l 


2)'(2n+-1)! 


, 


| —4(1—»)(1 — 27) Tomson oF 21 am+2n} 


l 
(2m—3)'(2n+1)!' 


‘ 


2n ‘ ss é 3 \ 
"Hi, (3—4n) Tom +2n 2— 21 am +2n-1 y J xm +n-2} 


2m—1 


where 


x 


| us du i use~ 2" 


: (s J, = a du 
sinh 2u—2u J sinh 2u—2u 


| u* 


isinh 2u 


{ ute? 
J,= | 


‘ {sinh 2u—2u 





‘ (4.28) 

The integrals J’, J, have been tabulated by Howland and Stevenson 
(11) and the asymptotic values shown to be the same as those for J,, J, 
respectively. The coefficients (4.28) are the same as (4.16) except that /,, 
J, are replaced by J’, J’, and the signs of J, are changed. The iterative 
procedure is now completely determined and by using (3.9), (4.15), (4.27) 
the series (2.3) may be constructed as far as required. 

It would appear that the expansions (4.11), (4.12), together with the 
integrals /,, J,, I’, J‘, are fundamental to the theory of thick plates. Their 
use would have avoided the approximation for the factor |2u+-sinh 2u]-! 
introduced by Sneddon ((6) 267) in order to evaluate integrals of the form 


u” sinh 
————_______ zu J,(au) du. 
2u+sinh 2u cosh 
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It is easily shown that 


“ © 


u” : | + , 
| —_.——— sinhzuJ,(wu) du > ete smth _ pv t2m+l Pr nei 

2u+sinh 2u a, (2v+2m+ 1)! 
0 h 


. x 


u” > > Tay +t v+2m Pv 

| saTeinn ae oeh term) de pa (vt omy Eee 

0 m= 

5. The prescribed stresses on the faces. Significance of the di- 
vergent terms 


As has already been stated, in the construction of D, an extra term must 
be added to account for the prescribed stresses across z = -+-h when these 
differ from zero. Suppose that the given boundary conditions on z = +h 
are first separated into their two parts corresponding to symmetric and 
antisymmetric systems, i.e. suppose 

% = f(w)+g(a), 
to = +j(w)+Kw) (5.1) 
on z= +h. 

The part of D, accounting for the symmetric stresses 22 = f(w), 
ta = +j(m~) on z = +h may be found without difficulty in terms of the 
Hankel transforms of f(w) and j(a) following the above procedure. It is 
of interest, however, to look more closely at the part of D, accounting for 
the antisymmetric stresses 2 = + 9(w), 2 = k(w) on z = +h, as this 
again gives rise to divergent terms in the integrands of the stress functions. 
Let . “ 
j= ( wy(a)J,(Ear) da, k= ( wk(o)J,(Ea) da. (5.2) 

0 0 
Provided modifications are made to ensure convergence, the required 
stress functions are then 


x 


= | EAcoshéa(éw) dé, fb = | EDsinh ézJ,(ar) dé, (5.3) 
0 0 
where A and D satisfy (4.19) with their right-hand sides replaced by £@, 
tk respectively, i.e. 


ain Ae (G cosh £h-+-k sinh €h}, 


. (g{ (1 —2n)cosh &h—€h sinh fh]+k 2(1— n)sinh £h —£h cosh €h]}}. 


= 2S 


D= 


(5.4) 
Expanding J,(éa), J,(€a) in powers of &, integrating the right-hand sides 
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of (5.2) term by term, and pasate in sili it is found that 


x 


£A cosh £zJ,(Ea) - | wg(a) da + 


sae eterdet + OF), 


Jr 


eR 
(5.5) 
ED sinh ézJ,(Ear) = Be. =e | wg(ar) dar + O(€ lp 
263;3 | 1 
0 
‘ » > 3p 
| sare | wre) dew"? + O10) (5.6) 
making use of (4.11), (4.12) ; 

The first terms on the right-hand sides of (5.5), (5.6) represent rigid 
body displacements and may be removed from their respective integrands 
in (5.3) without affecting the stresses. As before, the second terms in each 
of the right-hand sides of (5.5), (5.6) do not affect the stresses across 

th and so, after multiplication by e-*" to secure convergence at the 
upper limit, they may be removed from their respective integrands in 
(5.3). The integrals (5.3) are then convergent and determine the required 
part of D,. 

It will be noticed that the divergent terms, apart from those represent- 
ing rigid body displacements, given in (5.5), (5.6) cancel those in (4.22), 
(4.23) (for r = 0) if 


Zz 


( wq(u) da — 2h*p}'(1—7). 


0 


Also the resultant force P over the two faces is 
x « 
2 ( 2roy(w) dak = 4x [ wy(a) dak. 
0 0 
Hence, the integrals determining the stress functions for D, contain no 
divergent terms and therefore would not have needed any modification if 


_p 
My) = deceie (5.7) 
Sar( 1 — n) 

Now @_, corresponds to an isolated force at the origin equal to 87(1— »)k. 
Therefore the condition (5.7) is satisfied when the resultant force on the 
spherical cavity is balanced by the resultant force on the faces. More- 
over, in all succeeding terms D, (r > 1), p\*” = 0 since from (3.10) 


“nr 


Pi = Pi = py = 0. 

However, if there is a force resultant over the spherical boundary which 
is not balanced by the tractions on the faces, it must be balanced ‘at 
infinity’. For completeness, therefore, it is desirable to investigate the 
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behaviour of the stress functions (4.24), (4.25) for r = 0 as a + «©. These 
stress functions may be split up as follows 

d ™ ba + dy 
where 
d,, \fh + 3(1 Pe [7 Jy(Ear) —e-"* dé 4- 


4 


of _ b+ Wp, 


(1) 

+3(1— ne at Jo(&ar) — 1+ }é%ar2e-€"} dé, 
0 

‘ l 

lEA cosh &2 - 


+3()— p12 _ a. ) Pu. ‘| 
z(t 3(1—) a 3(1 1) gaol) dE, 
(1—ny(l—2 | | etsieor—« én dé - 


—3(1—y)(1 apt : {i a { Jy(Ear) — 1 + 1 £2 aye Eh dé, 
‘ aise 
|e) sinh €: -|ahe—(1- ~9)(1—2) = 1 | af. 


+3(1—n)(1— 29) Bh v7 \lem dé. 


Now Dougall has shown ((1) 136) that 


x 


P l 
| gialée)—e#) dé = -log(! mn) 
0 


and 


| gateolsor)—1-+ beterte 4 dg = fot 


Hoe(si) 


—_ ~__ 
[Thus it may be shown that the stresses wz, wom corresponding to the 
combined solution @, +p, are 


-_ 


ot = 3(1 h 


(1)y 
3(1— \-+n) Pt tog (F) —301 — Pe 'z 


p{tdzs 
+(g—f)- er. : 


_ 
wm 


h 
stant between z = this [ 27m a@2dz 
; —h 
the force on the cavity 
5092 .50 


Therefore, the resultant force in the z-direction over a nike ™ = con- 


—82(1—n)p\Ph?, which balances 
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However, the flexural couple 
h 
ie; 
2h 


a 


asaw—> ©. 


The term —(1+7)(1—yn)pPhlogm is the same as that given by the 
solution for an isolated force at the origin of magnitude 87(1—7y) in 
the approximate theory of thin plates; this gives some confirmation to 
the above analysis. 

Following the methods of Dougall (1) it may be shown that, as wa > 0, 
the stress functions ¢,, ys, merely annul the effect of the term h?p{!/r in 
the stress function ¢ determining D,. 

Hence, to sum up, p\!’ + 0 only in problems in which there is a force 
resultant Pk over the spherical boundary. If this force is balanced by 
the tractions on the faces the divergent terms in the integrands determining 
D, are annulled. If, on the other hand, either the whole or part of Pk is 
balanced at infinity the divergent terms must be removed according to 
the foregoing procedure and, as has been shown, an additional flexural 
couple must then be applied. In both cases p{"*) = 0 (r i 2.2...) e00 
so none of the integrals in the stress functions determining the succeeding 
terms D,, (ry -- 1) contain divergent terms (apart from those merely repre- 
senting rigid body translations). 


6. Convergence 
Sufficient conditions for the convergence of the symmetric part of D,,., ,, 
1.e. : 


2n4+1f »,(2r+1) (2r+1), ’ 
> h (Pon Pon t+dan Ponsa} 
n 1 


outside a sphere of radius p, < p are 
pa?) <M, igi) < ae (6.1) 


for some constant k, and where A, = p,/h < Ay. It will be assumed for 
the present that these are in fact satisfied. Also the asymptotic values 
of the coefficients given in (4.16) are 


(2m+2n-+-1)! 


9m 1. 9n)\! 
(en a Ss (2m-+-2n)! 
(2n+-1)!(2m— 1)! 22m+2n 


am ™ 75 D\'/9m 1)! Dams2n-1’ 
(2m— 2)!(2n4-1)!2 


2n+l hp 





2m 


(2m+ 2n-+-2)! ons] —(2m+-2n+1)! 
. ae - _ a ~ a 
(2m — 1)! (2n-+ 2)! Q2m+2n+1 am" (2m — 2)! (2n+ 2)! 22m+2n 


(6.2) 


9 ’ 
2n G, ~ 
om 





for large values of n or m. 
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Use will be made of the identity 


< 
= (v+r—1)! 


Tir)! (iz| <1), (6.3) 


(1—z)-” 
r-0 

where vy is a positive integer. 

Now for sufficiently large values of n, the moduli of the coefficients 
given in (4.15) are each less than their asymptotic values (6.2) multiplied 
by some constant (a say) > I. 

Therefore, from (4.15), 
‘. ~ {  (2m+2n-+ 1)! Az ‘ (2m + 2n)!A?"—! | 


](2r +2) ; ao = smindiina erent: siibe senineesingummeiadaeiaiente -} 
re = \(2n+- 1)! (2m— 1)! 22m +2" * (2m — 2)! (2n-4- 1)! Q2m+2n-1) 


. (2n+2) a,\-*- 2 aaa 
kody ( 5 by (6.3) 
— Akad, (2n-+4 2) 
° (2 A,)?" +3 
Similarly, it may be shown that 
4had,(2n + 3) 


(2r+2)) . —_—____——__.---- . 
Mons} (2 A, )®*+4 


Hence the series 


x 
, ; 2nfpi2r+ 2) 1 (2r + 2) ' 
Ds, .2 : 2h mln +1 Pansat Mean+i Ponsa} 
n 


for the symmetric part of D,,,,, is convergent inside a sphere of radius 
at least (2—A,)h which is at least greater than h, since A, < A < 1. Also, 
when the above inequalities for )*{?)) and |m3¥/?'| are substituted in 
(3.9), it is seen that, for sufficiently large values of n, 

pr <A, | < EAM, (6.4 
where k’ is some constant and A, < AZ < Ay. 

Now (6.4) are of the same form as (6.1) except for the factor n? which 
leaves the ensuing proof unaffected. Hence, provided the series defining 
D, are convergent, those defining D,,,, (r = 1, 2,3,...) are also convergent 
throughout the whole region occupied by the material. The series defining 
D,, (r = 1,2,3,...) are convergent within a sphere of radius at least half 
the thickness of the plate, and for points outside the sphere the values 
of the stress functions may be deduced from their integral forms. 

The convergence of the series used in dealing with the antisymmetric 
displacements may be proved similarly. 

For simple conditions of stress at infinity and stress-free faces, the series 
defining D, will usually be finite, in which case (6.1) would certainly be 
satisfied for r = 0. Hence, all the series used would converge at points 
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sufficiently near the cavity. In other cases, the convergence of D, would 
have to be examined separately. 


As stated in section 2, no consideration is given to the convergence of (2.3). 


7. Example (i). The plate in tension 

In this example the faces and cavity are considered to be stress-free 
and a uniform tension, parallel to the faces, applied at infinity, i.e. 

R.=0 om r=p, R,=0 on z= +h, 
oo>T as wD. 
It is easily verified that D, p, +m, satisfies the conditions at in- 
finity and also gives R, = 0 if 
mi”) 2. 

1+ 7 

In the numerical work, the value of 7 was taken to be } and the diameter of 
the cavity to be half the thickness of the plate, i.e. A, = 4. The coefficients 
(3.10), (4.16) were calculated to five significant figures and are tabulated in 
(9). By means of (3.9), (4.16) the terms of the series (2.3) were constructed as 
far as D, and the peripheral stresses yy and 60 around r = p were calculated 
at 15° intervals of 6. Comparison may be made with the approximate 
solution, obtained by neglecting the effect of the faces, that is D,+D,. 

The values of yy differ from those given by the approximate solution by 
about 7 per cent, at most. The values of #9 are somewhat closer to the 
approximate values; 99 varies around the cavity ina manner similar to that 
of the peripheral stress in the analogous two dimensional problem (see (10)). 


8. Example (ii). The plate in flexure 
In this example, the faces and cavity are again considered to be stress- 
free and a uniform flexural couple M applied at infinity, i.e. 


R 0 om r=~p, 2 th, 


r 


h 


= =e 
| wo dz 
2h. 


zum dz > M 


a as w > OD. 

It is easily verified that D, = h-4/{p,+mPp,} satisfies the conditions 
at infinity if sa 3M ities: —(1—2n)M 
7 2hA(1- n)’ . 2h(1+-) 
and also R, = 0. 
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The construction of successive pairs of stress functions defining the 
terms D, (r = 1, 2,3) of the series (2.3) follows as in the previous example 
except that formulae (4.27) must be used in place of (4.15) since the stress 
system is now antisymmetric. 

Comparison is made, as before, with the approximate solution. It is 
found that the latter solution is more accurate than in the previous 
example; the error being everywhere less than 3 per cent. for both 60 and 
yy around the cavity. Again the behaviour of 0 is similar to that of the 
peripheral stress in the analogous two-dimensional problem (11). 


9. Conclusion 

The analysis has been developed to a point at which it is possible to 
find the solution to any axially symmetric, torsion-free problem in the 
given region. In the two worked examples, the deviation from the infinite 
medium solutions are remarkably small as compared with the correspond- 
ing two-dimensional problems. This result might be expected for the 
following reason: when A, = 1, in the two-dimensional problems the strip 
is completely severed; in the three-dimensional problems, however, there 
are merely singularities on the faces. Thus, in the latter case, there would 
be a much less violent variation in the stresses as A, increases from zero 
(i.e. in the infinite medium solution) than in the former. 
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Values of the peripheral stresses around the cavity 
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Ay = 0 represents the approximate solution in which the effect of the faces is 








TORSION-FREE STRESS SYSTEMS 


REFERENCES 
. J. DouGay, Trans. Roy. Soc. Edin. 41 (1904) 129. 
ibid. 49 (1914) 895. 
A. E. Green, Phil. Trans. A 240 (1948) 561. 
and T. J. Witumore, Proc. Roy. Soc. A 193 (1948) 229. 
ibid. 195 (1948) 533. 
. I. N. Sneppon, Proc. Camb. Phil. Soc. 42 (1946) 260. 
7. Cutn-Bine Line, Quart. App. Math. 10 (1952) 149. 
ibid. 13 (1955) 381. 
. N. Fox, Ph.D. Thesis (London, 1958). 
. R.C. J. HOwLanpn, Phil. Trans. A 229 (1930), 49. 
and A. C. STEVENSON, ibid. 232 (1933) 155. 
. M. J. Boussinese, Application des Potentials (Paris, 1885). 
. B. GALERKIN, Compt. Rend. 190 (1930) 1047. 
P. F. Papkovitrcnu, ibid. 195 (1932) 513, 754. 
. E. STERNBERG, R. A. EuBanks, and M. A. Sapowsky, Proc. of the First U.S, 
National Congress of Applied Mechanics (Chicago, 1951). 
I. N. SNepvon, Fourier Transforms (McGraw-Hill, 1951). 
. E. T. Wurrraker, Modern Analysis (Cambridge, 1902). 
R. N. KaurMan, Soviet J. App. Math. and Mech. (P.M.M.) 22 (1958) 451. 





HEAT TRANSFER BY LAMINAR FLOW FROM A 
ROTATING SPHERICAL CAP AT LARGE 
PRANDTL NUMBERS 


By C. B. BAXTER and D. R. DAVIES (University of Sheffield) 
[Received 4 June 1959] 


SUMMARY 

An approximate method, described recently by Davies, for calculating heat 
transfer from a rotating disk is extended to the problem of calculating the distribu- 
tion of rate of heat transfer by laminar flow from a heated, rotating, spherical cap 
for an arbitrary, meridional distribution of surface temperature and for large Prandtl 
numbers. By comparing numerical results, for the case of constant surface tempera- 
ture, with those previously obtained for a heated rotating disk of equal area, an 
approximate evaluation is made of the effect of curvature of a rotating surface of 
revolution on heat transfer. 


1. Tue technique, used by Davies (1), for calculating heat transfer from 
a rotating disk is applied to the case of a heated, rotating spherical cap. 
The axes of spin and symmetry are assumed to coincide, and the speed 
of rotation is assumed to be sufficiently rapid to permit application of 
the boundary layer approximation; frictional heating is neglected. 


The surface temperature of the cap is assumed to be independent of 
position along circles of latitude and an arbitrary function of position 
along meridian ares. Using spherical polar coordinates, we denote distance 
from the centre O of the spherical surface to a point P by r, the angle 
between OP and the axis of spin by 6, the velocity component at P along 
the meridian by u, the stream function by #, the thermal diffusivity by 
x, and the temperature by 7. The temperature equation in von Mises 


form is then 1 eT a eT 
ma oe = Ka (ru |, (1) 
sin?6 06 op cys 
The exact distribution of w is not available, but Howarth (2) has applied 
a Karman-Pohlhausen momentum integral method to construct an ap- 
proximate solution. Following his analysis closely, we assume the poly- 


nomial form u(aQsin 8 cos 6X2)? = f(Z), (2) 

where {(Z) = 0-2357Z—0-5Z?+-0-5Z3— 0-4428Z4+0-2071Z°; a is the 

radius of the spherical surface, A is a function of @, and Z = 2/5, where z 

denotes normal distance from the surface, and 6 is the Karman boundary 

layer thickness, given by (v/Q)!A(@), v being the kinematic coefficient of 

viscosity and Q the angular velocity. This is similar to the polynomial 
(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 2, 1960] 
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form used by Howarth. It satisfies all the conditions specified by Howarth, 
and leads to a closer agreement with the exact values of skin friction, 
calculated by Cochran (3), when @ is small (3 per cent deviation). The 
polynomial is not, of course, valid in the vicinity of the equatorial plane. 
However, the heat is conveyed along the meridian in the direction of @ 
increasing, and we are therefore able to treat the problem of heat flow 
from a spherical cap, extending from @ = 0 to @ = 80° (approximately), 
which is the limit of the calculation made by Howarth. 

The associated expression for the stream function is easily seen to be 

| a®(Qv)? sin?@ cos 8A3|-! = | f(Z) dZ, (3) 
0 
and the dependence of u on % calculated numerically, from (2) and (3). 
The analytic relationship between u and %, required in order to obtain 
solutions of (1), is in general complex. However, the Lighthill (4) form 
of approximation is available in the asymptotic case of very large Prandt! 
numbers. In this case the first terms on the right-hand sides of (2) and 
(3) dominate, and we find 
u(aQ sin 8 cos 6 A?)-! = 0-69y4{a?(Qv)! sin?6 cos 6 A3}-+. (4) 
This form is found to be a very good approximation to the actual relation- 
ship between u and & up to Z = 0-08, the maximum deviation being 5 
per cent. 
If we now substitute (4) into (1) and make the transformation 
8 
d ( (sin @)2(cos 6)4At dé, 
0 
the basic equation becomes 
eT —. (vs a): 
ch Ob\" Ob 
where ¢c = 0-69Kv-!Qia3. 

2. Equation (6) is now identical in form with the basic temperature 
equation discussed previously by Davies (1, equation 9) and a solution, 
giving the distribution of local rate of heat transfer, Q, obtained for any 
prescribed meridional distribution of surface temperature. If 7,(¢,) de- 
notes the temperature of the surface at 4 = ¢», corresponding to 0 = 4p, 
and 7, that of the surrounding stationary fluid, we find that the distribu- 
tion of Q(¢_) is given by 
Q(bo) = (0°69)*(1-5)'2-! sin(27 3) P(3)kKQ4y-4e4(sin 0)(cos A,)# x 
$e 
} o-#)-TGo)—To] 44. (7) 


“di 


d 
dd, 
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In the case of constant surface temperature, the total rate of heat 
transfer, Z, from a spherical cap, 0 < 6 < @,, is given by 


E\ k(T, —T,)Qhv-ta?|-1 = 2-6811(3)sin(27/3)d} 04, (8) 


where & is the thermal conductivity of the fluid, and ¢, is the value of ¢, 
given by (5), corresponding to 6,. In the asymptotic case of very large o 
and very small @, we find that the value of E given by (8) is only 6 per 
cent greater than the value given by the exact solution for a disk (1). 


3. An approach can now be made to the problem of evaluating approxi- 


mately the influence of curvature on the rate of heat transfer from a 
rotating surface of revolution, by comparing the result obtained for a 
rotating spherical cap FE, with that obtained for a rotating disk of equal 
area E,,. We easily find that 


Ey Ey = $4[(2/38)Ag0 —cos 6)|-, (9) 


where A, is the value of A at @ = 0, and E,/E,, is taken to be 1-0 at @ = 0. 
The numerical values of this ratio are shown in the table for 0 < 6 < 80°. 


TABLE | 


; rate of heat transfer from a rotating spherical cap Ex 
Calculated values of - f fer f a ps 


rate of heat transfer from a rotating disk of equal area E,, 





0 ° 20 | 40 | 60 





EslEp 1°00 0°98 0°94 | 0°89 o'79 


Although the greater curvature of the spherical surface may be expected 
to lead to a higher rate of heat transfer, the radial and meridional velocity 
components in the cases of the disk and cap, respectively, clearly play 
dominant roles. If we calculate the ratio of the meridional component 
of velocity for a cap at a given Z, given by Howarth, to the corresponding 
radial component of velocity for a disk, given by Cochran, we find that 
as @ increases this ratio decreases progressively from 1-0 at @ = 0. This 
accounts qualitatively for the corresponding decreasing dependence of 
E,/E, on 6. If, of course, the curvature of a spinning surface of revolution 
is intermediate between that of a sphere and a flat surface, the values of 
Ex E, will fall between 1-0 and those given in the table. 

Finally, it should be emphasized that, since an exact calculation of the 
distribution of the meridional velocity u is not yet available, our calcula- 
tion of heat transfer is based on the approximate evaluation of the 
distribution of u given by Howarth, and is consequently of a similar 
approximate nature. 
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STRESS IN A DYNAMICALLY LOADED 
HELICAL SPRING 


By N. J. DURANT (25 Borough Road, Bridlington, E. Yorks.) 
| Received 23 October 1958 ; revise received 12 March 1959} 


SUMMARY 
In this paper a general expression is derived for the stress in a dynamically loaded 
helical spring. The solution has been effected by the Laplace transform and any 
time dependent force may be substituted for the force applied to the spring. 
A numerical example is worked for a spring of a particular mechanism—the swash 
plate pump. 


1. Introduction 


THE particular spring chosen for investigation is one of a number in a 


swash plate pump, a mechanism which as far as the analysis is concerned 


is shown in Fig. 1. 











lend 
Fic. 1. R, rotating shaft. s, swash plate integral with R. B, ball 
bearings in a circular race. RP, rocking plate supported by pistons 
and springs. P, piston, spring, and telescoping shrouds. 


The springs of such a pump are pre-loaded in compression for four 
reasons: (i) for the return motion of the pistons; (ii) to keep the pistons in 
contact with the rocking plate; (iii) to keep the rocking plate in contact 
with the ball or roller bearings, and (iv) to maintain contact of the bearings 
with the swash plate. 


[Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 2, 1960] 
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The springs are shrouded so that should a spring fracture, no additional 
damage to the mechanism would result. The springs therefore support 
the pistons, rocking plate, shrouds, and bearings. These masses supported 
by the springs and oscillating with them are considered in the analysis. 

As each spring is subjected to the same motion and to the same force, 
it is sufficient to consider the behaviour of one spring. The angular 
velocity of the rotating shaft varies from about four thousand to twelve 
thousand revolutions per minute according to the design of the pump and 
the amount of fluid to be pumped, but usually a pump is designed to 
operate at one maximum number of revolutions. 

The material of the spring wire is considered, as usual, to be homogeneous 
and isotropic. The angle which a coil makes with a generator, and the 
radius of a coil are taken to be constants, neither assumption being strictly 
true. It is known, however, that the stresses resulting from changes which 
must occur are small, if the angle of the helix is small. As a result of pre- 
loading and because the springs are confined within the limits of their 
pre-loaded length they are always in compression. The analysis shows 
that the rotational speed of the shaft has to be large for the occurrence 
of resonance. 


2. Solution 
We have to solve the one-dimensional wave equation 
$ 4 a © Weoes x - 0), (1) 
ox” c* od 
in which € is the axial dynamical displacement of an element of the spring, 
« is measured from that end of the spring remote from the applied force, 
c is the velocity of the stress waves in the spring wire, / is the maximum 
fitted length (pre-loaded length) of the spring, and ¢t is the time measured 
from a particular epoch. 
This eigenvalue problem is most easily solved by means of the Laplace 
transform. The subsidiary equation corresponding to (1) is then 


“(ye = 0, 
dx? \c 


where é is the Laplace transform of £; the solution of this equation is 
— — A cosh px/c+ Bsinh pz/c, (3) 
where A and B are arbitrary constants. 
From the initial condition, € = 0 when x = 0, and hence A = 0. The 
second initial condition is 


- 
= 


2 


fw) _ _ Fw, (4) 


t OX) p=} 
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where M is the total mass affecting the spring, k the spring constant, and 
F(t) is the dynamical force applied to the spring. 
The subsidiary equation corresponding to (4) is 
dé | R 
S| = —Fip). 


6 
dx\,_; 


Mp%—kl 


From (3) and (5) we obtain a 
. F(p)sinh pac 
(kB. pl.c)}B cosh pl/c—(pl'c)sinh plc} 
where B = ki? Me?. 
If F(t) = bks(1 cos wt) 
, ; a }h:sw? 
its transform is F(p) = —* nil , 
p(p* + w*) 
where s is the travel of the piston and w is the circular frequency. 
Substituting (8) in (6) and using the Inversion Theorem, the expression 
for the dynamical displacement is 
, Y +4 x 5 
sw*B e“ sinh Ax ¢ 
doi J (A*L c)(A2+-w?)| 8 cosh Al ce — (Al c)sinh A/c} 
y-—ix 
The integrand is a single-valued function of A with singularities (all 
poles) at 


dx. (9) 


A = @; A= +t: A= titca,/l (n 
where x, are the zeros of «tan « B. 
The residues at the poles are quoted below. If the residue at a singu- 
larity is Ry, then 


Ry_ 9 = 1/(Bw?l), 


c COS wt SIN wx /C 





~ lw*® B cos wl ¢ +-(wl /c)sin wc’ 


R of? S cos «, ctl sin a, x1 
nm atom “\e x,, (a2 — wl? c?)(a2 +B? B)cos a,, 


n=l 





By Cauchy’s theorem the value of the integral in (9) is 277 times the sum 
of the residues at the poles within the contour, therefore 
{x Be COS wt sin wax /c 
& = le ~~ = —_— 
- Lawl Beosal/c¢+(wl/c)sin wl (ec 


—2a(2) S cos x, ct/l sin a, 2/1 , (10) 


as | x, (a2 — wl? /c*)(a2 + 8? —B)cos a,, 








which is the expression for the dynamical displacement for 0 < x < l, 
t> 0. 
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If (10) is differentiated once with respect to x we obtain the expression 
for the dynamical strain, which if multiplied by kl gives the expression 
for the dynamical force, which is: 
o& 


Ox 


Lh { , 8 cos wt cos wx/¢ 
= ks se Rcteteeteeted SS 
- a 


8 cos wl/c+(wl/c)sin wl /c 


kl 


2 a ~] . o 
2a() See eremea er __| 
" 2 2]2 /¢2)( 2 2 
c a (a2 —w?l?/c?)(a? + 8? —B)cos a,,| 
Q@ say. (11) 
This completes the analytical solution. We now require the dynamical 
stress, to which must be added the static stress resulting from the pre- 


load in the spring. 


3. Numerical example 
We consider a spring and its oscillating masses having the following 
properties which will be used in the numerical example: 
Material Steel 
Modulus of elasticity E = 30.10° |lb/in? 
in shear ( = 12.108 |b in? 
Diameter of wire d = 0-08 in. 
coil = 0-66 in. 
Maximum fitted length of spring = 1-25 in. 
Travel of piston s = 0-50 in. 
Angle of helix = 7° 
Spring stiffness > = 25 Ib/in. 
Circular frequency of shaft = 4,200 r.p.m. = 1407 
radians sec. 
Weight of piston +, == 0-072 lb 
spring *» = 0-028 Ib 
| (rocking plate+shrouds+ bearings) w, = 0-070 1b 
Total oscillating weight 7 = (w,+w,/3+w,)t = 0-15 lb 
Mass density of steel = 0-283 /g 
Number of coils = 10 
Number of pistons ie | 
Pre-load in spring = 8lb 
c= ,(E/y) 202.103in. /see 
B = kl?/Mc* = 2-462.10-* 
wl 'c = 8664.10-7 


| To show the order of magnitude of the 
quantities involved. 


Also we have atana = —B ~ 0 hence a, = nz (n = 1,2, 3.,...). 


+ See Appendix. 
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From the order of magnitude of the quantities given above we can 
estimate the value of Q (equation (11)) without detailed numerical work. 
The coefficient of the third expression in (11) is a small quantity. The 
numerator cannot exceed unity, and the denominator is of the order 
(—)"an*, so that, for this particular case, this part of (11) may be neglected. 
To a close approximation we then have 


Q = $ks(1—cos wt cos wx/c). (12) 


The maximum value of this expression occurs when wt = 7 and x = 0, 
when Q = ks = 12} 1b. (13) 


This is very nearly the value of Q when wt = 7 and x = / because wl/c 


is small. 
The dynamical stress in the spring wire is given by 
o, = 8QD(1+sin «)/7d* 
= 46,040 lb in? 
and the statical stress resulting from the pre-load of 8 |b is 
op = 46,040 x 8/124 = 29,470 lb/in?®. 
The total stress which results from flexure and shear is 
o = 6,,+op = 75,510 Ib/in*. 
Equation (11) then shows that with the numerical values a, and wl/c 
resonance will not occur in this particular case. 


APPENDIXt 
In the list of properties we have taken one-third of the weight of the spring as 


its oscillating weight. We now investigate whether this is justifiable. 
If we consider the original homogeneous equation 


(A 1) 


then a solution of this equation satisfying the first boundary condition = 0 when 
’ dis ; 
slain € = asin qgxcos(qet +e). (A 2) 


The second condition: ay = 

Ot? | 2} 

will be satisfied if qitangl = m/M 

which, if we write gl = z, is ztanz = m/M. 

If m/M is small then, expanding tanz, we have 
24(1+ 427+...) = m/M 


+ I am indebted to the referees for suggesting writing the appendix. 
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and to a first order approximation we can write 


lm ) m 
"3M! M’ 
and thus 28 m/(M 4 4m). 


The period of the lowest frequency is 


2*( | 


Ty 27/cq 


(A7) 


Qari | {/M +4m 
=/{ m). 


c m 


Using the numerical data in section 3 we have M+ 4m = 0-150.9g,m = 0-028/9 hence 


T; 27r(le)./ (53571) 27r(1/c)2-3145. (A 8) 
If no approximation had been made in solving (3a) we would have 


ztan2 mM ‘ (AQ) 


The smallest root of this equation is z 


04328, hence from the formula 7, 2zl 'cz 
we find that 


T; 277 (L/c)2-3105. (A 10) 

The periodic times calculated either from (A 5) or from (A 9) are almost identical, 

so that the approximation of adding one-third of the spring mass to the other 

oscillating masses is justified in this particular case, The approximation of adding 
one-third of the spring mass would not be valid if the ratio m/)M were not small. 
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